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ABSTRACT  of  Final  Report  of  F49620 -01-1-0142 


Lattice  Boltzmann  equation  (LBE)  was  used  to  perform  direct  numerical  simulations 
(DNS)  and  large-eddy  simulations  (LES)  of  benchmark  problems  in  turbulence, 
scalar  mixing  and  reaction.  The  most  important  contribution  was  the  development 
of  LBE  theory  for  binary  mixing  which  can  be  extended  in  a  straight-forward 
manner  to  multi -scalar  mixing.  Computational  simulations  were  performed  to  verify 
that  the  desired  diffusive  effects  are  achieved  in  classical  mixing  problems.  For 
example,  it  was  demonstrated  that  the  shape  of  the  probability  density  functions 
during  binary  mixing  from  non-premixed  initial  conditions  were  precisely 
captured.  The  technique  for  handling  reactions  in  the  LBE  context  was  also 
demonstrated.  In  the  standard  one-dimensional  flame  propagation  problem,  the 
burning  rate  was  captured  accurately.  The  third  significant  contribution  was  the 
adaptation  of  multi -time- scale  relaxation  technique  to  LES.  Several  DNS  and  LES 
of  benchmark  turbulent  flows  (decaying  isotropic  and  homogeneous  shear,  square 
jet)  were  performed  to  establish  the  effectiveness  and  efficiency  of  LBE.  The 
decay  exponent  in  decaying  turbulence,  the  equilibrium  anisotropies  in 
homogeneous  turbulence,  and  the  spread  rates  in  square  jets  were  accurately 
calculated.  In  summary,  a  solid  foundation  for  future  LBE-LES  calculations  of 
turbulent  combustion  has  been  established. 
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Chapter  1 


General  Introduction 


Physically  accurate  and  computationally  viable  calculation  procedures  for  turbulent  com¬ 
bustion  are  required  for  a  variety  of  propulsion  applications  ranging  from  improving  combus¬ 
tor  performance  to  estimating  the  concentrations  of  environmentally-harmful  trace  chemi¬ 
cals  and  strategically-crucial  signature  producing  species  precisely.  Despite  the  vast  strides 
in  computational  capabilities,  direct  numerical  simulation  (DNS)  of  turbulent  combustion 
for  device-size  computations  remains  well  beyond  the  capabilities  of  current  computers, 
as  well  as  those  for  the  foreseeable  future.  It  is  now  a  widely  held  view  that  large-eddy 
simulation  (LES),  in  conjunction  with  the  probability  density  function  (PDF)  method  for 
the  subgrid  scales,  offers  the  best  hope  for  computing  practical  turbulent  combustion  prob¬ 
lems.  The  moment  methods  are  generally  unsuited  for  reacting  flows  because  of  the  strong 
non-linearity  introduced  by  the  reaction-rate  term.  The  LES-PDF  technique  is  an  attrac¬ 
tive  choice,  as  it  combines  precise  simulation  of  the  problem-dependent  large-scale  motion 
with  the  reasonably  accurate  PDF  method  for  the  universal  small-scale  fluctuations.  While 
much  less  computationally  intensive  than  direct  numerical  simulations  (DNS),  a  clean  LES 
simulation  with  the  spectral  cut-off  well  in  the  inertial  range  (as  dictated  by  theory)  still 
can  be  very  expensive.  Many  applications  of  interest  to  the  Air  Force  involve  multi-phase 
flows  (sprays  and  atomization),  preferential  diffusion  requiring  the  computation  of  multi- 
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component  diffusivity,  and  complex  geometry.  For  these  problems,  the  LES  method  can  be 
very  tedious,  to  the  point  of  being  computationally  intractable.  In  the  PDF  methodology, 
advection  and  reaction  terms  appear  in  closed  form  but  diffusion  processes  require  closure 
modeling.  Physically  accurate  diffusion  models  tend  to  be  computationally  expensive,  and 
the  simpler  models  can  be  quite  inaccurate.  As  currently  practiced,  the  evaluation  of  the 
reaction-rate  term  can  take  up  to  80%  of  the  computational  effort  per  time  step.  Reduced 
chemistry  models  are  used  typically  to  alleviate  this  computational  burden.  The  reduced 
chemical  kinetics,  while  adequately  capturing  the  slow  reactions  and  major  species,  may  be 
inaccurate  for  trace  species  that  could  be  environmentally  harmful  or  signature  producing. 
Most  traditional  reduction  techniques  assume  that  many  trace  species  are  in  partial  equilib¬ 
rium,  leading  to  errors  in  estimation  of  the  concentration  of  these  species.  Even  in  the  more 
sophisticated  reduction  techniques,  such  as  the  Intrinsic  Low-Dimensional  Manifold  (ILDM) 
method,  the  fast  eigen-species  are  assumed  to  be  in  equilibrium,  leading  to  less  accurate  es¬ 
timation  of  these  species  than  the  slow-evolving  major  species.  While  the  LES-PDF  method 
is  the  best  technique  available  currently,  there  is  a  lot  of  room  for  improvement. 

In  this  report,  we  describe  the  early  development  of  a  new  methodology  for  calculating 
turbulent  combustion  based  on  the  more  fundamental  Boltzmann  equation  rather  than 
Navier-Stokes  equation.  The  ultimate  goal  of  the  proposed  work  is,  in  effect,  a  large  eddy 
simulation  of  the  Boltzmann  equation  in  which  the  (molecular)  velocity  distribution  function 
is  discretized  in  phase  space  with  a  lattice.  More  succinctly,  we  call  our  approach  the  large 
eddy  simulation  (LES)  of  the  lattice  Boltzmann  equation  (LBE).  This  LBE-LES  technique 
has  several  physical  and  computational  advantages  over  the  LES-PDF  approach. 

1.1  Lattice  Boltzmann  vs.  Navier-Stokes  equations 

The  Boltzmann  equation  is  potentially  a  better  hydrodynamic  platform  for  LES  calculations 
of  turbulent  combustion  than  is  the  Navier-Stokes  equation.  The  advantages  lie  both  in 
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improved  physical  accuracy  and  better  computational  characteristics.  The  physical  accuracy 
aspects  will  be  discussed  in  detail  later. 

Computational  advantages.  The  computational  philosophies  of  the  two  methods  are 
vastly  different-  While  the  LBE  method  leads  to  very  simple,  explicit  and  frequent  calcula¬ 
tions  at  each  node,  the  NS-based  calculation  at  each  node  is  complex  (often  implicit)  and 
occurs  less  frequently.  At  first  glance,  it  would  appear  that  the  advantage  achieved  by  LBE 
due  to  the  simplicity  of  the  calculation  would  be  offset  by  the  higher  frequency  of  compu¬ 
tations,  rendering  the  two  methods  equal.  That  is  not  to  be.  What  sets  the  two  methods 
apart  is  the  locality  of  the  calculations  at  each  node.  The  LBE  scheme  leads  to  calculations 
that  are  not  only  simple  but  also  involve  only  local  interactions.  In  the  NS  methods,  global 
interactions  play  a  major  role,  resulting  in  a  severe  communication  penalty  and  a  marked 
reduction  of  computational  speed.  This  difference  becomes  even  more  significant  when  com¬ 
putations  are  performed  on  multi-processor  parallel  platforms.  The  NS-based  computations 
typically  scale  poorly  with  an  increase  in  the  number  of  processors.  The  LBE  computations, 
on  the  other  hand,  display  excellent  scalability  characteristics;  even  super-linear  scaling  has 
been  observed  in  some  very  special  cases.  The  main  advantages  of  the  LBE-based  methods 
are: 

1.  LBE  is  better  suited  for  large-scale,  especially,  parallel  computing. 

2.  LBE  is  ideally  suited  for  handling  multi-phase  flow  with  phase  transition  and  multi¬ 
species  mixtures,  where  multi-component  diffusivity  is  important.  Navier-Stokes 
solvers  can  be  computationally  too  expensive  for  these  flows. 

3.  LBE  can  handle  complex  geometry  with  relative  ease.  Even  computations  involv¬ 
ing  moving  boundaries,  e.g.,  rotating  turbine  blades,  can  be  handled  without  loss  of 
computational  speed. 
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In  our  own  work,  we  have  determined  that  LBM  is  a  factor  of  six  times  faster  than  a 
commercial  second-order  finite  difference  code  (FIRE)  for  isotropic  643  DNS  simulations. 
We  expect  this  factor  to  be  larger  for  larger-size  computations. 

1.2  Objectives  and  accomplishments. 

The  objective  of  the  research  is  to  develop  the  fundamentals  of  the  Lattice- Boltzmann  method 
to  enable  large-eddy  simulations  of  turbulent  combustion.  Our  ultimate  goal  is  to  compete 
with  NS  methods  in  their  strength  areas  and  provide  a  superior  computational  capability  in 
problems  that  are  not  easily  tractable  by  the  NS  method. 

The  areas  of  focus  of  the  completed  research  are: 

1.  Molecular  mixing,  especially  between  species  of  vastly  different  molecular  weights. 

2.  DNS  and  LES  of  turbulence. 

3.  Laminar  flames. 

4.  Chemical  kinetics  and  their  reduction. 

5.  Treatment  of  complex  geometry. 

During  this  period,  important  progress  in  the  areas  of  molecular  mixing,  laminar  flames, 
treatment  of  complex  geometry  in  laminar  flows,  DNS /LES  of  turbulence  in  simple  geometry 
and  chemical  kinetics  reduction  has  been  made.  In  the  area  of  molecular  mixing,  we  have 
developed,  from  kinetic  theory,  the  basic  LBM  model  for  multi-component  mixtures.  We 
also  demonstrated  that  the  LBM  mixing  model  can  reproduce  classical  results  and  obtain 
important  new  results  that  are  beyond  the  scope  of  continuum-based  methods.  LBM  can 
simulate  mixing  and  segregation  with  equal  facility.  In  the  area  of  laminar  flames,  we 
have  developed  computational  capability  that  accurately  calculates  the  burning  velocity 
in  a  simple  one-dimensional  premixed  flame.  Since  practical  combustors  are  typically  of 
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complex  geometry,  we  have  developed  preliminary  versions  of  complex  boundary  treatments. 
Reducing  chemical  kinetics  is  as  important  in  LBM  as  it  is  in  NS-based  computations.  For 
LBM,  we  have  evaluated  stochastic  chemistry  and  deterministic  chemistry  methods.  We 
have  concluded  that  stochastic  kinetics  models  are  not  yet  feasible  due  to  computational 
issues. 

An  incomplete  list  of  refereed  journal  articles  published  from  research  (partly  or  fully) 
supported  by  this  grant  is  given  below. 

1.  L.-S.  Luo  and  S.  S.  Girimaji,  Lattice  Boltzmann  model  for  binary  mixtures ,  Phys.  Rev. 

E  66(4):035301(R)  (2002).  [Addresses  Theory  of  mixing] 

2.  L.-S.  Luo  and  S.  S.  Girimaji,  Theory  of  the  lattice  Boltzmann  model:  Two-fluid  model 
for  binary  mixtures,  Phys.  Rev.  E  67:036302  (2003).  [Addresses  Theory  of  mixing] 

3.  H.  Yu,  L.-S.  Luo,  and  S.  S.  Girimaji,  Scalar  mixing  and  chemical  reaction  simulations 
using  lattice  Boltzmann  method,  Int.  J.  Computat.  Eng.  Sci.  3(l):73-87  (2002). 
[Addresses  random  mixing  and  chemical  reaction] 

4.  D.  d’Humieres,  I.  Ginzburg,  M.  Krafczyk,  P.  Lallemand,  and  L.-S.  Luo,  Multiple- 
relaxation-time  lattice  Boltzmann  models  in  three- dimensions ,  Philos.  Trans.  R.  Soc. 
London  A  360(1792)  :437-451  (2002).  [Addresses  fundamental  LBM  issues] 

5.  M.  Krafczyk,  J.  Tolke,  and  L.-S.  Luo,  Large-eddy  simulations  with  a  multiple-relaxation¬ 
time  LBE  model,  Int.  J.  Mod.  Phys.  C  17(l/2):33-39  (2003).  [Addresses  fundamental 
LES  issues] 

6.  P.  Lallemand  and  L.-S.  Luo,  Hybrid  finite- difference  thermal  lattice  Boltzmann  equa¬ 
tion,  Int.  J.  Mod.  Phys.  C  17(l/2):41-47  (2003).  [Address  heat-release  physics  model] 

7.  P.  Lallemand  and  L.-S.  Luo,  Lattice  Boltzmann  method  for  moving  boundaries,  J. 
Computat.  Phys.  184(2) :406-421  (2003).  [Addresses  complex  moving  boundary  is¬ 
sues] 
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Several  more  papers  are  undergoing  review  or  under  preparation. 


Chapter  2 


Main  Results  and  Inferences 


Turbulent  combustion  involves  three  important  physical  processes:  (i)  turbulent  stirring, 
(ii)  molecular  mixing,  and  (iii)  chemical  reaction.  The  first  two  constitute  the  process  of 
turbulent  mixing.  The  stirring  of  a  scalar  field  (temperature,  concentration,  etc.)  by  a 
turbulent  velocity  field  causes  gradient  steepening,  leading  to  enhancement  of  molecular 
mixing.  Since  turbulence  is  essentially  a  3-D  phenomena,  our  LBE  computation  of  tur¬ 
bulence  is  three-dimensional  in  nature.  Molecular  diffusion  process,  on  the  other  hand,  is 
essentially  one-dimensional  in  nature:  the  diffusion  is  along  the  direction  of  the  scalar  con¬ 
centration  gradient  which  is  orthogonal  to  the  iso-scalar  surface.  Hence,  the  LBE  studies 
of  molecular  mixing  have  been  conducted  in  three  and  two  dimensions.  Chemical  reaction 
is  a  local  (zero-physical  dimension)  process.  Hence  we  assess  the  suitability  of  the  LBE 
for  chemical  reaction  in  a  simple  1-D  flame.  In  our  work,  we  develop  and  test  the  LBE 
models  for  these  processes  individually.  Extensive  simulations  have  been  performed  leading 
to  important  results.  In  this  Chapter,  we  summarize  the  main  results  and  inferences. 

A  novel  lattice  Boltzmann  model  for  binary  mixing  is  proposed  and  discussed  in  Section 
2.1  (details  in  Appendix  B).  In  Section  2.2,  we  simplified  the  mixing  model  further  and  test 
that  model  in  non-premixed  binary  scalar  mixing  (details  in  Appendix  C  and  D).  In  Section 
2.3,  LBE  is  tested  for  one-dimensional  flame  propagation  through  a  premixed  mixture  of 
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propane  and  air  (details  in  Appendix  D).  Results  from  DNS  and  LES  of  turbulence  in 
isotropic  and  homogeneous  turbulence  are  presented  in  Section  2.4  and  2.5,  respectively 
(details  in  Appendix  E  and  F).  Finally,  a  LBE  model  for  axisymmetric  geometry  is  given  in 
Section  2.6  (details  in  Appendix  G).  Further  work  of  LBE-LES  on  investigation  of  mixing 
enhancement  using  square  jet  is  not  presented,  as  it  is  not  yet  complete. 

2.1  Lattice  Boltzmann  model  for  binary  mixtures 

2.1.1  Introduction 

The  present  day  lattice  Boltzmann  equation,  with  its  high  fidelity  physics  and  computation¬ 
ally  efficient  formulation,  is  a  viable  alternative  to  the  continuum  methods  for  simulating 
fluid  flows.  In  fact,  it  can  be  argued  that  for  many  complex  problems  involving  multi-fluid 
phenomena,  the  physics  can  be  captured  more  naturally  by  the  Boltzmann  equation-based 
methods  rather  than  Navier-Stokes  equation-based  methods.  There  is  a  growing  interest  in 
using  the  lattice  Boltzmann  equation  (1;  2;  3;  4;  5;  6;  7).  However,  most  existing  LBE  mod¬ 
els  for  multi-component  fluids  (8;  9;  10;  11;  12;  13;  14;  15;  16)  tend  to  be  somewhat  heuris¬ 
tic  and  make  the  single-fluid  assumption.  The  single-relaxation-time  or  Bhatnagar-Gross- 
Krook  (BGK)  approximation  (17)  is  used  in  most  existing  models  (10;  11;  12;  13;  14;  15;  16), 
restricting  applicability  to  unity  Prandtl  and  Schmidt  numbers. 

In  this  work,  we  develop  a  two-fluid  lattice  Boltzmann  model  that  is  based  on  kinetic 
theory  for  binary  mixtures.  Such  a  model  would  be  capable  of  (i)  simulating  arbitrary 
Schmidt  and  Prandtl  numbers,  and  (ii)  accurately  modeling  the  interaction  between  miscible 
and  immiscible  fluids.  We  follow  a  general  approach  within  the  framework  of  kinetic  theory 
for  developing  the  lattice  Boltzmann  models  for  multi-fluid  mixtures.  This  work  is  a  part 
of  our  ongoing  effort  to  set  the  lattice  Boltzmann  equation  on  a  more  rigorous  theoretical 
foundation  and  extend  its  use  to  more  complex  flows.  We  derive  a  discretized  version  of  the 
continuum  Boltzmann  equations  for  binary  mixtures.  The  extension  of  this  methodology 


2.1.  Lattice  Boltzmann  model  for  binary  mixtures 


13 


to  multi-fluid  mixtures  is  relatively  straight-forward. 


2.1.2  Kinetic  theory  of  gas  mixtures  and  lattice  Boltzmann  equation 

The  simultaneous  Boltzmann  equations  for  a  binary  system  are  (please  see  Appendix  B  for 
nomenclature) 


dtfA  +  £■ V/A  +  aA- Vf/A  =  Qaa  +  Qab,  (2.1a) 

dtfB  +  £■ V/B  +  Ob- Ve/B  =  Qba  +  Qbb,  (2.1b) 


where  QAB  =  QBA  is  the  collision  term  due  to  the  interaction  among  two  different  species 
A  and  B. 

With  the  BGK  approximation  (17;  18),  the  collision  integrals  Qcq  [a,  q  €  (A,  B)]  can 
be  approximated  by  the  following  linearized  collision  terms 


=  (2.2a) 

A(T 

j*  =  -^-[r  -  r?{0)],  (2.2b) 

where  and  /<TC(°)  are  Maxwellians 


o)  _ 


Tin 


0-(S-Uv)2/(2R<tT<t) 


(2nRcrT^D/2 

/«( 0)  _  n<T  c-(£-v*')2/(2R<rT*') 


(2n  R^)0'2 


(2.3a) 

(2.3b) 


where  D  is  the  spatial  dimension,  Ra  —  kg / ma  is  the  gas  constant  of  the  a  species,  kg  is 
the  Boltzmann  constant  and  ma  is  the  molecular  mass  of  the  a  species.  There  are  three 
adjustable  relaxation  parameters  in  the  collision  terms:  \a,  A?,  and  ACT<r  =  (nc/ncr)A?cr. 
The  first  Maxwellian  /a(°)  is  characterized  by  the  conserved  variables  of  each  individual 
species:  the  number  density  nc,  the  mass  velocity  ua,  and  the  temperature  Ta.  The  second 
Maxwellians  and  are  characterized  by  four  adjustable  parameters:  ua<;,  u?(T, 

TCTf,  and  TC(r. 


2.1.  Lattice  Boltzmann  model  for  binary  mixtures 


14 


For  the  isothermal  flows,  we  can  derive  the  lattice  Boltzmann  equation  by  discretizing 
the  model  equations  (2.1)  following  the  procedure  described  in  (4;  5)  (for  details,  please  refer 
to  Appendix  B): 

fa(xi  +  ea  5t,t  +  ^t)  —  fa(xi >*) 

=  J™  +  J?-FZ6t,  (2.4) 


where  the  self-collision  term  the  cross-collision  term  and  the  forcing  term  Fa  are 
given  by 


•C  =  “[.£- /£<0)] . 

■C  =  ~  “)■(“•  _  “<)> 

td  p  Kal 

■cut  ®a  ‘ 

Fa  ~  —Wap(j  r  t  , 


(2.5a) 

(2.5b) 

(2.5c) 


where  pa  and  pq.  and  ua  and  u?  are  the  mass  densities  and  flow  velocities  for  species  a  and 
?.  They  are  the  moments  of  the  distribution  functions: 


a >  =  ££  =  ££“”'  (2-6a) 

a.  a 

P<TU<7  =  faect  =  /a(0)eQ>  (2-6b) 

a  a 

and  p  and  u  are,  respectively,  the  mass  density  and  the  barycentric  velocity  of  the  mixture: 

P  =  pc  +  Pc,  (2.7a) 

pu  =  PpU<j  +  pcUc-  (2.7b) 


The  collision  terms  J ™  and  J£c  are  constructed  in  such  a  way  to  satisfy  the  mass  and 
momentum  conservation  laws. 

The  equilibrium  distribution  function  fa^  has  the  following  form  in  general: 


/q(°)  = 

/£(6q)  =  WaPa 


1  +  J^T^01  “  U^'  ^Ua  ~  U>> 


fa 

J  a 


a(eq) 


1  .  (eg  ~  y)  •  U  (eg  -u)2' 

RaT  2(RaT)2  ‘ 
where  coefficients  {u;Q}  depend  on  the  discrete  velocity  set  {eQ}. 


(2.8a) 

(2.8b) 
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2.1.3  Discussion  and  conclusions 

We  have  constructed  a  lattice  Boltzmann  model  for  binary  mixtures  with  several  important 
features.  All  the  modeling  issues  are  addressed  within  the  framework  of  extended  kinetic 
theory.  The  lattice  Boltzmann  model  then  is  derived  directly  from  the  continuous  kinetic 
model  equations  using  a  formal  discretization  procedure.  The  lattice  model  thus  inherits  the 
sound  physics  and  mathematical  rigor  incumbent  in  kinetic  theory.  This  rigor  is  in  contrast 
to  previous  lattice  Boltzmann  models  for  mixtures  (11;  12;  13;  14;  15;  16),  which  are  not 
based  directly  on  the  fundamental  physics  of  continuum  kinetic  equations.  These  models 
rely  on  fictitious  interactions  (11;  12)  or  heuristic  free  energies  (13;  14;  15;  16)  to  produce 
the  requisite  mixing.  (Many  defects  of  the  free-energy  models  (13;  14;  15;  16)  are  due  to  the 
incorrectly  defined  equilibria  (19).)  These  non-physical  effects  present  a  further  problem, 
since  they  are  not  easily  amenable  to  mathematical  analysis  (19;  20).  The  heuristic  elements 
of  the  previous  lattice  Boltzmann  models  (11;  12;  13;  14;  15;  16)  have  been  eliminated, 
resulting  in  a  physically  justifiable  model  that  is  simple  to  compute.  Further,  due  to  the 
close  connection  to  kinetic  theory,  the  derivation  of  the  hydrodynamic  equations  associated 
with  the  lattice  Boltzmann  model  is  simplified  significantly  and  rendered  mathematically 
more  rigorous.  The  derivation  of  the  hydrodynamic  equations  from  the  previous  lattice 
models  are  much  less  rigorous  (11;  12;  13;  14;  15;  16). 

The  second  important  feature  of  the  present  work  is  that  the  model  is  based  upon  a 
two-fluid  theory  of  binary  mixtures.  The  previous  models  (8;  9;  10;  11;  12;  13;  14;  15;  16), 
on  the  other  hand,  are  derived  from  a  simpler,  but  highly  restrictive,  one-fluid  theory.  In 
the  single-fluid  models  with  the  BGK  approximation  one  is  constrained  to  use  the  ad  hoc 
“equilibrium  velocity”  (11;  12) 

u(eq)  _  iTsP<ru(r  4~ 

( T<lPa  4"  TcrP<;) 

in  the  equilibrium  distribution  function  fa ^  in  order  to  satisfy  the  local  conservation  laws. 
As  a  result,  the  viscous  relaxation  process  and  the  diffusion  process  are  inseparable.  The 
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analysis  of  these  models,  therefore,  becomes  unnecessarily  tedious  and  cumbersome  (11;  12). 
The  models  with  free  energies  (13;  14;  15;  16)  do  not  yield  correct  hydrodynamic  equations 
(19;  20),  mostly  due  to  the  incorrectly  defined  equilibrium  distribution  functions  used  in 
these  models  (19).  Furthermore,  single-fluid  models  cannot  be  applied  to  mixtures  of  species 
with  vastly  different  properties.  In  the  present  two-fluid  model,  the  diffusion  behavior  is 
decoupled  from  viscous  relaxation.  The  diffusivity  is  determined  by  the  parameter  td  and 
the  physical  properties  of  the  mixture.  The  model  is  capable  of  simulating  either  miscible 
or  immiscible  fluids  by  changing  the  sign  of  ( td  —  1/2). 


2.2  Simple  model  for  binary  scalar  mixing 

While  the  binary  mixing  model  derived  in  the  previous  Section  is  general  and  mathemati¬ 
cally  rigorous,  it  may  be  difficult  to  implement.  In  this  Section,  we  will  develop  a  relatively 
simple  model  for  binary  mixing  by  adding  inter-particle  reaction.  This  model  is  amenable 
to  easier  implementation  and  testing.  (For  nomenclature  and  details,  please  see  Appendix 
C  and  D) 


2.2.1  Lattice  Boltzmann  equation  for  scalar  mixing 

Let  a  and  a'  denote  the  two  species  of  interest.  The  LBE  for  each  species  is 


n°(x  +  ea,t  +  1)  =  n%(x,t)  + 


(2.9) 


where  na  is  number  density  distribution  function  and  the  collision  operator  is  given  by 

i  jaa ' 

=  1  +  tSt  (2-io) 


mu 


includes  an  additional  term  J°a  that  reflects  interactions  between  the  two  species. 
We  use  the  following  number  density  equilibrium  distribution  function 


=  wana[l  4-  3(eQ  •  tf)  +  ^(ea  •  v!)2  -  ^ u ,2] 


(2.11) 
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with 


*-££<»*«•) 


The  binary  interaction  term  is  modeled  as 


(2.12) 


Jf  =  •  [V^  +  (s'  -  *0^1  =  -tf' 


(2.13) 


where  x°  and  ua  are  molar  and  mass  fractions  of  the  species  a 


x 


a  _ 


na  „  pa 


= 


n 


and 


n^q)  =  wan[  1  +  3(eQ  •  u)  +  ^(eQ  •  «)2  -  ^u2] 


(2.14) 


The  hydrodynamic  properties  of  the  above  equations  are  given  in  Appendix  C. 


2.2.2  Simulations  and  Conclusions 

We  performed  non-premixed  binary  scalar  mixing  computation  in  a  periodic  box.  This 
problem  epitomizes  the  molecular  mixing  issues  encountered  in  a  typical  non-premixed 
combustion  application.  Two  species  (presumably  fuel  and  oxidizer)  are  segregated  initially 
and  randomly  distributed  in  the  computational  domain,  which,  in  the  present  case,  is  a 
square  box.  The  two  species  are  labeled  generically  as  black  and  white.  A  typical  initial 
distribution  is  shown  in  Figure  2.1. 

The  macroscopic  velocity  is  set  everywhere  to  zero,  corresponding  to  a  pure  diffusion 
problem.  The  initial  values  for  the  number  densities  are  nb  =  1.0,  nw  =  0.0  in  region  of  the 
black  species  and  nw  =  1.0,  nb  =  0.0  in  region  of  white  species.  Assuming  homogeneity  of 
the  scalar  field,  periodic  boundary  conditions  are  used  in  all  directions. 

The  first  case  studied  is  the  mixing  of  two  fluids  of  equal  molecular  weights  and  number 
densities (mb  =  mw  =  1.0),  hence  of  equal  mass  density.  This  case  is  particularly  interesting, 
as  the  results  can  be  compared  directly  with  direct  numerical  simulation(DNS)  of  the  Navier- 
Stokes  equation  data  of  Eswaran  and  Pope  (21).  Figure  2.2(a)  shows  the  time  evolution  of 
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Figure  2.1:  Initial  number  density  distribution  for  both  equal  and  unequal  mass 
cases 


the  probability  density  function(pdf)  of  the  scalar  p: 


The  corresponding  DNS(21)  data  are  shown  in  Figure  2.2(b).  The  LBE  and  DNS  data  show 
excellent  qualitative  agreement. 

In  Figure  2.3,  the  time  evolution  of  the  root-mean-square(rms)  of  scalar  fluctuations 
(p1)  obtained  from  LBE  is  compared  with  that  from  DNS(21).  They  agree  quite  well.  We 
also  performed  the  simulation  for  an  unequal  mass  case(mh  =  2.0,  mw  =  1.0).  The  detailed 
results  are  given  in  Appendix  D. 

In  conclusion,  we  present  a  binary  mixing  model  and  use  it  to  simulate  scalar  mixing.  In 
the  case  of  equal-density  species  mixing,  well-known  results  from  a  continuum  Navier-Stokes 
simulation  are  reproduced.  The  true  advantage  of  the  LBM  can  be  seen  from  the  mixing 
simulations  of  species  of  different  molecular  weights.  The  results  appear  quite  encouraging. 
Such  simulations  are  very  difficult  with  continuum  based  methods. 
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2.3  Lattice  Boltzmann  equations  for  reacting  flow 

2.3.1  Introduction 

We  simulate  simple  one-dimensional  flame  propagation  through  premixed  propane  and  air. 
The  problem  studied  is  identical  to  that  of  Yamamoto(22),  but  the  physical  combustion 
field  is  sightly  different.  The  simplifying  assumptions  invoked  in  this  study  are: 

•  There  are  no  external  forces  in  the  field. 

•  The  flow  is  incompressible. 

•  The  chemical  reaction  (heat  release)  does  not  affect  the  flow  field;  the  temperature 
and  concentration  fields  are  decoupled  and  solved  separately. 

•  Nitrogen  is  inert. 

•  The  transport  properties  axe  constants  (not  functions  of  temperature). 

•  Viscous  energy  dissipation  and  radiative  heat  losses  are  negligible. 

•  A  simple  one-step  reaction  is  considered 

C3H8  +  502  =►  3 C02  +  AH20 

and  the  over'" ill  reaction  rate  is  given  by 

U>ov  =  Kov  CczHs  C()2  e  E/KT 

where  Cc3Hs ,  Cfo2>  kov  and  E  are  the  concentrations  of  propane  and  oxygen,  the 
reaction  coefficient,  and  the  effective  activation  energy,  respectively. 

The  details  of  modeling  this  problem  are  given  in  Appendix  D.  Here  we  will  only  present 


numerical  results. 
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2.3.2  Simulations  and  conclusions 

A  schematic  of  the  simulated  flow  is  shown  in  Figure  2.4.  For  this  simple  case,  the  back¬ 
ground  flow  maintains  both  the  pressure  and  velocity  fields  uniform  in  space  and  time.  A 
heat  source  is  placed  at  a  location  close  to  the  inlet  to  ignite  the  mixture.  Once  ignition 
is  achieved,  the  heat  source  is  removed.  At  subsequent  time  the  flame  propagates  to  the 
right. 

In  Figure  2.5  the  flame  position  is  shown  as  a  function  of  time.  The  flame  location 
is  identified  as  the  position  with  the  highest  reaction  rate  at  any  given  time.  The  linear 
variation  of  flame  location  with  time  (in  Figure  2.5)  indicates  that  the  flame  propagates  at 
a  nearly  constant  rate.  The  burning  velocity  is  Sl  =  0.11  m/s,  which  compares  extremely 
well  with  the  value  obtained  from  experiments  for  a  propane-air  flame(23). 

Figure  2.6  shows  the  reaction  rate  profile  in  the  reaction  zone  as  time  evolves.  Simu¬ 
lations  indicate  that  flame  behavior  is  sensitive  to  the  magnitude  of  the  heat  source.  The 
premixed  reacting  flow  simulations  produce  results  that  are  in  good  agreement  with  known 
data.  Based  on  these  simulations,  we  believe  that  LBM  can  perform  adequately  for  more 
complex  chemical  reaction  calculations. 


L - H 


Figure  2.4:  A  schematic  illustration  of  a  simple  1-D  reacting  flow 
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2.4  DNS  and  LES  of  decaying  homogenous  isotropic  turbu¬ 
lence  (HIT) 

2.4.1  Introduction 

We  now  examine  the  ability  of  LBM  to  simulate  turbulence  cascade  which  is  responsi¬ 
ble  for  stirring  a  scalar  field.  We  perform  three  type  of  simulations:  faVLBE-DNS  of 
decaying  HIT  in  inertial  and  rotating  frame  of  reference.  The  decay  exponents  for 
the  kinetic  energy  k  and  the  dissipation  rate  e  are  computed  and  compared  with  correspond¬ 
ing  NS-DNS  results.  The  low  wave-number  scalings  of  the  energy  spectrum  are  studied.  The 
effect  of  rotation  on  the  kinetic  energy  decay  is  investigated,  (b)  LBE-LES  of  decaying 
HIT  inertial  frame  of  reference.  We  compute  the  kinetic  energy  decay,  energy  spec¬ 
trum,  and  flow  structures  using  LBE-LES.  By  comparing  LBE-LES  results  with  the  corre¬ 
sponding  LBE-DNS  results,  we  observe  that  LBE-LES  accurately  captures  large  scale  flow 
behavior.  We  find  that  the  optimal  Smagorinsky  constant  value  for  LBE-LES  is  smaller 
than  the  traditional  value  used  in  NS-LES  approaches,  (c)  LBE-LES  vs.  NS-LES.  We 
carry  out  a  comparative  study  of  the  LBE-LES  and  NS-LES  of  decaying  HIT.  We  show  that 
the  LBE-LES  simulations  preserve  flow  structures  more  accurately  than  the  NS-LES  coun¬ 
terpart.  This  accuracy  is  due  to  the  fact  that  some  history/non- local  effects  are  inherent  in 
the  LBE  subgrid  closure. 

2.4.2  Simulation  results 

LBE-DNS  of  decaying  isotropic  turbulence 

The  variation  of  decay  exponent  (n)  with  Taylor  Reynolds  number  (  Re*)  in  various  sim¬ 
ulations  is  shown  in  Fig.  2.7.  The  dependence  of  n  on  Re*  obtained  by  the  LBE-DNS  is 
very  similar  to  that  observed  in  NS-DNS  calculations  (24).  The  values  of  n  obtained  in  the 
present  work  agree  well  with  the  experimental  and  numerical  NS-DNS  data.  In  Figure  2.8, 
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we  can  see  that  the  low  wave-number  spectra  scale  as  E{k)  ~  «4  or  E{k)  ~  k2  depending 
upon  initial  spectrum.  This  dependence  of  low-wave  number  scaling  on  the  initial  spectrum 
is  in  exact  agreement  with  the  results  reported  in  (24;  25). 


Rex 


Figure  2.7_  Dependence  of  the  decay  exponent  n  =  1  /{C2  -  1)  on  initial  conditions  and 
ReA.  The  quantity  C2  is  depicted  m  the  figure  in  stead  of  n.  Solid  lines  represent  NS-DNS 
data  from  (24),  and  symbols  correspond  to  the  LBE-DNS  results  of  the  present  work.  For 
the  128  resolution,  •:  «rms  =  0.0064,  kmin  =  1,  kmax  =  8,  and  *  =  0.01  (r  =  0.53);  A: 
Urms  -  0.021,  fcmin  _  8,  kmax  =  16,  and  u  «  0.00167  (r  =  0.505);  o-  uTms  =  0  022  k  ■  =1 

*max  =  8,  and  I,  «  0.00167  (r  =  0.505).  For  the  643  resolution  (x):  uims  =  0.022  k'l  =  4 
and  kmax  —  8,  and  v « 0.00167  (r  =  0.505).  •  «ns  ’  min 


LBE-LES  of  decaying  isotropic  turbulence 

Our  first  effort  is  to  determine  the  appropriate  Smagorinsky  constant  ( Csm )  for  LBE-LES. 
Figure  2.9  shows  that  the  kinetic  energy  decays  with  different  Smagorinsky  constant  values 

for  both  323  and  643  cases.  We  find  that  Csm  =  0.1  is  the  best  choice,  and  we  shall  use  it 
in  all  calculations. 

In  Fig.  2.10,  we  compare  the  instantaneous  flow  structure  of  uz(i,  j,  k  =  N/2,  t')  ob¬ 
tained  by  the  LBE-LES  with  that  by  LBE-DNS.  As  shown  in  Fig.  2.10,  the  LBE-LES 

appears  to  capture  the  flow-field  structure  quite  adequately  even  with  a  coarse  resolution 
of  323  grid  points. 
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Figure  2.8:  Energy  spectra  for  two  cases  of  1283  above  with  initial  energy  concentrat¬ 
ing  in  the  range  of  1  <  «  <  8  at  t'  =  0  and  t'  =  0.022  respectively.  The  ini¬ 
tial  energy  spectra  are  given  by  E(k)  =  0.038k4 exp(-0.14K2)(upper  plot)  and  EU)  = 
0.038k  exp(— 0.14  k2) (bottom  plot)  respectively. 
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Figure  2.9:  Kinetic  energy  decay  with  different  Csm  and  resolution.  Solid  line  for  LBE- 
DNS  (128  )  and  symbols  for  LBE-LES  (o:  643,  Csm  =  0.1;  A:  643,  Csm  =  0.17;  x-  323 
Csm  =  0.1;+:  32s,  Csm  =  0.17). 

LBE-LES  vs.  NS-LES 

We  compare  the  LBE-LES  with  the  NS-LES  results  at  323  resolution.  The  difference  can 
be  seen  from  contours  of  the  instantaneous  flow  field  uz(i,  j,  k  =  JV/2,  t')  at  three  different 
times,  as  shown  in  Fig.  2.11.  By  Comparing  the  LBE-LES  results  (left  column)  and  the 
NS-LES  results  (right  column)  with  the  LBE-DNS  results  (center  column),  we  observe  that 
LBE-LES  preserves  the  flow  structure  better  than  the  corresponding  NS-LES.  The  323 
LBE-DNS  contours  shown  here  are  obtained  by  truncating  the  1283  LBE-DNS  data. 

2.4.3  Summary  and  conclusions 

First,  we  conduct  the  direct  numerical  simulations  by  using  the  LBM.  The  well  known 
power-law  decay  of  the  kinetic  energy  is  reproduced.  The  decay  exponents  obtained  in  the 
LBE  simulations  are  in  good  agreement  with  the  results  from  experimental  measurements 
and  NS-DNS  calculations.  Second,  we  conduct  a  comparative  study  of  the  LBE-LES  and 
the  LBE-DNS.  Comparisons  of  643  and  323  LBE-LES  with  1283  LBE-DNS  show  that  the 
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Figure  2.10:  Contours  of  the  instantaneous  flow  field  uz(i ,  k  =  N/ 2,  £').  LBE-DNS  vs. 
LBE-LES  with  different  resolutions.  The  323and  643  LBE-DNS  contours  shown  here  are 
obtained  by  truncating  the  1283  LBE-DNS  data 
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LBM-LES,  t'=0.006  LBM-DNS,  t' =0.006  NS-LES,  t’=0.006 


Figure  2.11:  Contours  of  the  instantaneous  flow  field  uz(i,  j,  k  =  N/ 2,  //) .  The  LBE-LES 
and  NS-LES  with  a  resolution  of  32s  compared  to  the  LBE-DNS  with  a  resolution  of  1283 
at  three  different  times.  The  323  LBE-DNS  contours  shown  here  are  obtained  by  truncating 
the  1283  LBE-DNS  data 
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large  scale  motion  is  captured  well  by  LBE-LES.  Finally,  we  compare  both  the  LBE-LES 
and  NS-LES  results  with  the  corresponding  DNS  results.  Our  comparisons  indicate  that 
LBE-LES  has  better  capability  to  preserve  flow-field  structure  than  NS-LES.  This  work 
further  establishes  the  LBE  method  as  a  viable  computational  tool  for  turbulence  stirring 
simulations. 

2.5  Study  of  homogenous  turbulence  shear 

2.5.1  Introduction 

We  will  first  verify  the  applicability  of  the  lattice  Boltzmann  method  (LBM)  in  the  DNS 
of  steady  homogenous  turbulence  and  then  focus  on  the  new  physics  in  flows  subjected  to 
time-dependent  shear.  For  details,  please  refer  to  Appendix  F.  Time-dependent  shear  can 
simulate  the  effect  of  mixing  in  wakes  and  mixing  layers  where  flow  is  stirred  by  vortices. 

2.5.2  Homogenous  shear  flow 

The  computational  domain  is  shown  in  Fig.  2.12.  An  uniform  1283  mesh  is  considered. 
Figure  2.13  shows  the  development  of  the  normal  components  of  anisotropy  tensor  (fry). 
LBM  agrees  with  the  experimental  data  and  numerical  results  very  well  (26;  27).  Fig.  2.14 
compares  the  present  612  (diagonal  component)  with  the  result  obtained  by  Jacobitz  et  al. 
(28).  The  agreement  is  again  very  good.  From  the  above  comparisons,  we  conclude  that 
DNS  using  LBM  can  recover  classical  results  accurately. 

2.5.3  Homogenous  turbulence  subjected  to  time- dependent  shear 

Here  we  let  shear  vary  with  time  and  investigate  the  effects  of  the  variation  on  the  devel¬ 
opment  of  turbulence.  The  periodic  variation  is  meant  to  simulate  the  effect  of  periodic 
stirring. 

Figures  2.15  and  2.16  show  the  evolution  of  turbulent  kinetic  energy  for  the  low  ( ^=0.5) 


2.5.  Study  of  homogenous  turbulence  shear 


30 


Figure  2.12:  Computational  domain. 


Figure  2.13:  The  results  of  anisotropy  tensor  obtained  from  DNS  using  LBM,  Navier-Stokes 
equations,  and  experiment.  St  is  non-  dimensional  time 
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Figure  2.14:  The  results  of  anisotropy  612  obtained  from  DNS  using  LBM  and  Navier-Stokes 
Equations  by  Jacobitz  et  al.  (28)  with  the  initial  values  of  Re>  =  44.72  and  Sk/e  =  2.0.  St 
is  non-  dimensional  time 

and  high  (|  =1.0)  frequency  cases  (  S  is  the  shear  rate  and  cj  is  the  frequency  of  the  variation 
of  shear).  When  frequency  is  low,  the  turbulence  kinetic  energy  k  increases  after  an  initial 
decay  while  for  the  high  frequency  case,  k  decreases.  Thus,  there  exists  a  critical  frequency 
at  which  turbulence  goes  from  growth-mode  to  decay-mode.  The  knowledge  of  this  physics 
is  important  for  evaluating  mixing  enhancement  strategies.  More  details  of  the  analysis  of 
this  phenomenon  are  given  in  Appendix  F. 

The  details  of  the  development  of  turbulence  anisotropy  tensor  are  given  in  Figs.  2.17 
and  2.18.  In  the  low  frequency  case,  the  anisotropy  in  turbulence  is  fully  developed.  The 
mean  values  of  6n  and  622  are  about  0.15  and  -0.1  respectively.  In  the  high  frequency  case, 
the  mean  values  of  bn  and  622  are  about  zero.  Thus,  turbulence  has  no  opportunity  to 
develop  as  in  the  high  frequency  case. 

Figures  2.19  and  2.20  show  the  evolution  of  production  over  dissipation  ratios.  There 
exists  negative  production  during  some  period  of  evolution.  For  the  low  frequency  case,  the 
average  value  of  P/e  is  above  one,  and  for  the  high  frequency  case,  this  value  is  around 
zero.  This  explains  why  the  turbulent  kinetic  energy  evolves  differently  in  these  cases. 
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Figure  2.15:  Evolution  of  the  turbulence  kinetic  energy  k  for  ^  =  0.5  case,  ko  is  the 
turbulence  kinetic  energy  at  the  initial  stage;  Smaxt  is  non-dimensional  time. 


Figure  2.16:  Evolution  of  k  for  ^  =  1.0  case 
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Figure  2.19:  Evolution  of  production  over  dissipation  ratio  j  for  ^  =  0.5  case 


Figure  2.20:  Evolution  of  production  over  dissipation  ratio  j  for  ^  =  1.0  case 
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2.5.4  Comparison:  DNS  vs.  RANS  models 

We  compare  DNS  results  with  the  Reynolds  averaged  Navier-Stokes  (RANS)  turbulence 
models.  We  use  3  different  models:  LRR-IP,  LRR-QI,  and  SSG  models  (29).  In  the  high 
frequency  case,  DNS  and  the  models  predict  decay  of  k,  as  shown  in  Fig.  2.21.  However 
the  models  predict  much  slower  decay  rates.  In  the  low  frequency  case,  the  models  fail  to 
predict  the  growth  of  k ,  as  shown  in  Fig.  2.22. 


0  20  s  (40  60 


Figure  2.21:  DNS  vs.  RANS  models  :  kinetic  energy  evolution  for  ^  =  1.0  case 

2.5.5  Conclusions 

DNS  of  homogenous  turbulent  shear  flow  (constant  shear)  is  performed  using  the  lattice 
Boltzmann  method.  DNS-LBM  results  agree  well  with  DNS-NS  and  experimental  data.  In 
DNS  of  homogenous  turbulence  subjected  to  periodically  varied  shear,  a  critical  frequency, 
which  determines  whether  k  decays  or  not,  is  found.  We  observe  the  negative  production 
of  k ,  which  means  that  there  is  a  mechanism  to  take  kinetic  energy  away  from  turbulence. 
The  comparisons  of  the  results  of  DNS  and  turbulence  models  show  that  the  models  fail  to 
predict  much  of  the  observed  behavior. 
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Figure  2.22:  DNS  vs.  RANS  models  :  kinetic  energy  evolution  for  j  =  0.5  case 

2.6  The  LBM  for  axisymmetric  flows  in  cylindrical  coordi¬ 
nate 

In  this  Section,  we  extend  the  LBM  to  a  cylindrical  coordinate  system  and  present  lat¬ 
tice  Boltzmann  models  for  axisymmetric  flows  for  both  the  vector(momentum)  field  and 
scalar(temperature,  concentration,  etc.)  field.  This  is  accomplished  by  inclusion  of  a  source¬ 
like  term  or  a  body  force-like  term(30)  in  the  lattice  Boltzmann  equations. 

By  introducing  a  body  force  term  ga( r,  t)  in  the  lattice  Boltzmann  equation  with  single¬ 
relaxation  collision  operator  —  the  so-called  BGK  model  (17) —  we  obtain 

fair  +  ea6t,t  +  St)  =  fa(r,t)  -  i[/a  -  f^]  +  5tga(r,t),  (2.15) 

For  the  cylindrical  system,  we  find  that  body  force  term  ga  can  be  expressed  as 

3po  1  c2  dp'  pc,  , 

ga  =  Was  +  -^waea-aL+-  —  —  +  Z-u-Vur  +  waea-aL  (2.16) 
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,  ds  ,  ds 

a‘  =  s,TTz'  a’  =  kT& 


2.6.1  Simulation  and  conclusions 

We  verify  our  model  using  a  pipe  Poiseuille  flow  as  it  amenable  to  an  exact  solution.  In 
Fig.  2.23,  the  flatness  corresponding  to  Eo / /i2  ~  constant  demonstrates  that  this  model  is 
second  order  in  space.  In  the  discussion,  E2  is  L2  norm  error  and  h  is  grid  size. 
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Figure  2.23:  Errors  between  numerical  and  analytical  solutions  computed  by  formula  G.48. 
G  =  2.89351e  -  05,  r  =  1.0 


In  conclusion,  we  present  a  LBM  model  for  axisymmetric  flows.  The  scheme  has  been 


proved  to  be  second  order  accurate  in  space. 
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Chapter  3 


Conclusions 


The  Lattice  Boltzmann  method  (LBM)  is  an  innovative  computational  tool  for  calculating 
complex  fluid  flow.  This  method  is  based  on  the  more  fundamental  Boltzmann  equation 
rather  than  the  Navier-Stokes  equation.  Due  to  its  firm  founding  in  kinetic  theory,  LBM 
offers  several  key  physical  advantages  for  computing  chemically  reacting  flow.  Most  im¬ 
portantly,  it  offers  the  possibility  of  modeling  the  processes  of  advection,  mass  transfer, 
and  scalar  and  thermal  transport  in  one  single  unified  description.  Furthermore,  chemical- 
reaction  computation  at  the  sub-continuum  level  is  easily  amenable  to  LBM  because  of 
the  its  kinetic  nature.  These  attributes  are  significant  advantages  over  continuum-based 
methods,  where  various  processes  are  modeled  with  individual,  often  incompatible,  models. 
Progress  made  in  the  research  under  this  grant  clearly  demonstrates  the  potential  of  this 
approach;  however,  much  more  needs  to  be  done  before  this  method  can  be  used  routinely 
for  large  scale  turbulent  combustion  calculations.  In  future  research,  we  shall  continue  to 
perform  fundamental  research  necessary  to  realize  that  goal. 

The  research  was  performed  principally  at  Texas  A  &  M  University  with  Prof.  Sharath 
S.  Girimaji  as  the  principal  investigator  (PI).  The  National  Institute  of  Aerospace  (NIA) 
Hampton,  Virginia,  was  a  sub-contractor,  with  Dr.  Li-Shi  Luo  as  the  co-PI.  We  also  col¬ 
laborated  with  the  leading  research  groups  in  the  world,  as  necessary,  to  advance  LBM  as 
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a  viable  computational  tool  for  turbulent  combustion. 


Chapter  4 


Appendix 


Appendix  A 


Background  of  Lattice  Boltzmann 
Method 


Nomenclature 

a<n> 

Coefficients  of 

a1' 

Acceleration  due  to  external  force 

Speed  of  sound 

&Ct 

Discretized  particle  velocity 

f 

Distribution  function 

/(°) 

Equilibrium  distribution  function 

hW 

Hermite  polynomial 

L 

Linearization  operator 

Q 

Collision  operator 

Pa/? 

Second  moments  of  / 

R 

Gas  constant 

T 

Temperature 

U 

Fluid  velocity 

P 

Fluid  density 

£ 

Particle  velocity 

/ 

After  collision  state 

A 

Relaxation  time 

T 

Normalized  relaxation  time 

Wa 

Weight  function 

Historically,  Lattice  Boltzmann  Equation  was  developed  from  lattice  gas  automata  (1; 
2;  3).  However,  the  lattice  Boltzmann  equation  can  be  directly  derived  from  the  continuous 
Boltzmann  equation  (4;  5;  6;  7),  which  can  be  written  in  a  general  form  for  a  iV-component 
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mixture  (8): 


dtf  +  e-Vf  +  ai-V^fi 


r  f *  (i + w  *)(i + *3  f '*)]*/*& 

N 

Ew./’i. 

3= 1 


(A.l) 


where  i9j  =  —1  0,  +1,  for  Fermi-Dirac,  Maxwell-Boltzmann,  and  Einstein-Bose  statistics; 
dtp  =  27 rcr  sin  x^Xj  the  quantities  a  and  x  are  the  differential  scattering  cross-section  and  the 
angle  of  deflection,  respectively;  p  =  f(x,  £*,  t)  and  f*  =  f(x,  £'\  t)  are  the  distributions 
before  and  after  collision  for  the  ?-th  species,  and  a 1  is  the  acceleration  due  to  external 
force  fields  or  a  mean-field  interaction.  Derivation  of  LBM  equations  from  the  Boltzmann 
equation  involves  three  steps  that  are  briefly  described  here. 

Linearization.  For  continuum  systems  or  even  neutron  scattering  processes,  it  is  suf¬ 
ficient  to  use  the  linearized  collision  term  L(g )  to  approximate  the  nonlinear  collision  term 
Q(f,  /)  as  /  =  /(°)(1  +  g)  and  only  keep  linear  terms  in  g,  where  /(°)  is  the  equilibrium 
solution,  i.e.,  fW)  =  0.  Function  g  can  be  expanded  in  some  appropriate  orthogonal 

basis,  such  as  the  generalized  Hermite  polynomials  {H^}  (9), 


9  = 

71=0 

where  £0  :=  (£  -  «),  a(0)  =  1,  Oo}  =  0,  and  a^j  =  (Pap/p  -  8ap), 


(A.2) 


Pa3  —  j  ^OaCopf  d£, 


The  linearized  collision  term  can  be  written  in  terms  of  {H^}.  For  Maxwell  molecules, 
the  linearized  collision  operator  becomes 

OO 

L(g)  =  J]  Ana(")  •  H<B\  L(H<B>)  =  AnH<n).  (A.3) 

71=0 

Reduction.  By  approximating  the  spectrum  {An}  of  L  by  a  finite  number  of  distin- 


46 


guished  eigenvalues,  the  linearized  collision  operator  becomes 
N  oo  N 

LN(g)  =  '£xn^-H^-ujN  £  aW.HW^fAn  +  ^aW.HW-^.  (A.4) 

n—0  n—N+ 1  n=0 

By  setting  N  =  0,  the  collision  term  is  reduced  to  the  simplest  linear  collision  model,  the 
celebrated  Bhatanaga-Gross-Krook  (BGK)  model: 

4  [/-/<»>],  A  =  i. 

Discretization.  We  shall  only  briefly  describe  the  derivation  of  a  lattice  Boltzmann 
model  from  a  continuous  kinetic  equation  (4;  5;  6;  7).  For  the  sake  of  simplicity  without 
losing  generality,  we  use  the  BGK  equation  for  a  single-component  system  without  the 
forcing  term  and  rewrite  it  as  the  following 

A/  =  ~[/-/(0)],  Dt~dt  +  t-V.  (A.5) 

Integration  along  characteristics  followed  by  Taylor  series  expansion  to  first  order  in  time 
leads  to 

f(x  +  Z,  t  +  St)~  /(*,  Z,  t)  =  I  [f(x,  £,  t)  -  /(°)  (x,  Z,  i)]  ,  r  :=  A  (A.6) 

Further  discretization  in  phase  space,  followed  by  appropriate  approximation  of  the  equi¬ 
librium  distribution  and  low  Mach  number  assumption  leads  to  the  current  popular  form 
of  LBM.  For  a  simple  single  species  problem  mixing,  this  leads  to  the  following  lattice 
Boltzmann  equations: 

fa(x  +  ea5t,t  +  5t )  -  fa(x,t)  =  (/a  -  /£0))  ,  (A.7) 

where  r  is  the  normalized  relaxation  time  which  determines  the  viscosity  and  fip  is  the 
discretized  equilibrium  distribution  function: 

(ea  -u)-u  (ea  ■  u )2 

RT  2(RT)2  ' 

(A.8) 


The  fluid  density  and  flow  momentum  are  moments  of  the  distribution  given  by 


p=J2f°‘  =  J2  ^0)>  PU  =  Y1  =  IZ  e<*/£0)>  (A-9) 

a  a  a  a 

In  the  above  equations  {wa}  are  uniquely  determined  by  the  discrete  velocity  set  {eQ}.  It 
is  important  to  point  out  that  since  the  LBM  is  intrinsically  multi-dimensional,  there  is  no 
essential  difference  between  two-  and  three-dimensional  models  -  the  only  difference  is  the 
number  of  discrete  velocities. 

The  most  popular  models  of  the  LBM  for  one-component  are  D2Q9  and  D3Q19  models 
for  2D  and  3D,  respectively,  as  shown  in  Figure  A.l  and  A.2. 


Figure  A.l:  D2Q9  model:  9  discrete  velocities  on  a  square  lattice 


Then  the  equilibrium  distribution  function  for  isothermal  fluids  is  given  as 
fi0)  =  wap[  1  +  3-^(eQ  ■  u)  +  i(eQ  •  u )2  -  \u2] 


(A.10) 


in  which,  for  D2Q9  model,  the  discrete  particle  velocities  ea  and  the  weighting  factor 
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Figure  A.2:  D3Q19  model:  19  velocities  lattice  in  3D 


wa(a  =  0, 1, 2,  •  •  •  ,  8)  are 


^Ql  —  \ 


(0,0), 

(cos[(o!  —  l)7r/2],sin[(a  —  l)7r/2]), 


a  =  0 

a  —  1,2, 3,4 


(A.11) 


(cos[(a  -  5)7t/2  +  7t/4],  sin[(a  -  5)7r/2  +  7r/4]),  a  =  5, 6, 7, 8 


and 


4/9,  a  =  0 

wn  =  <  1/9,  a  =  1,2, 3,4  (A.12) 

1/36,  a  =  5, 6, 7, 8 

respectively,  and  for  D3Q19  model, the  discrete  particle  velocities  ea  and  the  weighting 
factor  wa(a  =  0, 1, 2,  •  •  •  ,  18)  are 


[  (0,0),  a  =  0 

<  (±1,  0,  0)c,  (0,  ±1,  0)c,  (0,  0,  ±l)c,  a  =  1-6 

k  (±1,  ±1,  0)c,  (±1,  0,  ±l)c,  (±1,  ±1,  0)c,  a  =  7-18. 


(A.13) 


L 
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and 

t 

1/3,  a  =  0 

wa  =  \  1/18,  a  =  1,2,...,  6  (A. 14) 

1/36,  a  =  7, 6,...,  18 

s. 

respectively. 

The  sound  speed  is  cs  =  1/ >/3(5x/<5t)  with  5X  being  the  lattice  constant  of  the  underlying 
square  lattice.  The  macroscopic  quantities  such  as  mass  density  p  and  mass  velocity  u  are 
given  by 

P  =  £/a  (A.15) 

a 

pU  =  ^Tfaea  (A.16) 

a 

Further  points  of  note  are:  (i)  By  means  of  the  Chapman-Enskog  analysis  (multiple 
scale  expansion),  we  can  derive  the  hydrodynamic  equations  for  single-species  fluid,  (ii) 
This  LBM  formulation  can  be  extended  for  multi-species  systems  as  shown  in  (7).  (iii) 
The  lattice  Boltzmann  equation  can  be  constructed  as  a  fully  discrete  dynamical  system 
without  referring  to  the  Boltzmann  equation,  but  relying  to  some  first  principles  of  the 
kinetic  theory.  The  LBM  in  this  setting  is  called  the  generalized  LBE  method  or  multiple- 
relaxation-time  (MRT)  LBE  (10;  11;  12;  13).  Clearly  the  MRT-LBE  is  the  fully  discrete 
truncated  counterpart  of  the  linearized  Boltzmann  equation. 
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Lattice  Boltzmann  Model  for 
Binary  Mixtures 


a 

Nomenclature 

Acceleration  due  to  external  force  c 

Peculiar  velocity  of  molecule 

Das 

Binary-diffusion  coefficient 

Daa 

Self-diffusion  coefficient 

&OL 

Discretized  particle  velocity 

J 

Linearized  collision  operator 

/ 

Distribution  function 

/( 0) 

Equilibrium  distribution  function 

k 

Thermal  conductivity 

ks 

Boltzmann  constant 

m 

Mass  of  molecule 

n 

Number  density 

Q 

Collision  operator 

p 

Pressure 

R 

Gas  constant 

Sij 

Strain  rate 

T 

Temperature 

0 

Solid  angle 

u 

Fluid  velocity 

p 

Fluid  density 

4> 

Mass  concentration 

V 

Molar  concentration 

* 

Particle  velocity 

1 

After  collision  state 

A 

Relaxation  time 

T 

Normalized  relaxation  time 

wa 

Weight  function 

V 

Viscosity 

e 

Internal  energy 

St 

Time  step 

P 

Stress  tenser 

Q 

Heat  flux 

PD,PT,PM,v'M 

Relaxation  parameters 
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B.l  Introduction 

In  many  practical  flows  involving  pollutant  dispersion,  chemical  processing,  and  combustor 
mixing/reaction,  mass  and  momentum  transport  within  multispecies  fluids  plays  an  impor¬ 
tant  role.  For  these  flows,  it  can  be  difficult  to  construct  continuum-based  models  from  first 
principles.  Further,  these  flows  typically  involve  complex  geometry  and/or  multiple  phases 
making  computation  with  continuum-based  models  quite  complicated.  Therefore,  for  these 
flows,  there  is  a  growing  interest  in  using  the  lattice  Boltzmann  equation  (1;  2;  3;  4;  5;  6;  7). 

In  general,  for  computing  fluid  flow  of  any  type,  the  lattice  Boltzmann  equation  (LBE) 
(1;  2;  3;  4;  5;  6;  7)  offers  several  computational  advantages  over  continuum-based  methods 
while  retaining  the  flow  physics  intact.  Although  the  origins  of  the  modern  lattice  Boltz¬ 
mann  equation  (LBE)  can  be  traced  back  to  lattice-gas  automata  (8;  9;  10),  the  new  LBE 
models  are  free  of  some  well-known  defects  associated  with  their  predecessors.  Recent  works 
have  unequivocally  established  that  the  lattice  Boltzmann  equation  is  in  fact  connected  to 
kinetic  theory  (4;  5;  6)  and  completely  consistent  with  the  fundamental  conservation  prin¬ 
ciples  governing  fluid  flow  (11;  12;  13).  In  these  Appendix,  a  priori  derivation  of  the  lattice 
Boltzmann  equation  from  the  parent  continuous  Boltzmann  equation  is  developed  (4;  5;  6). 
The  Navier-Stokes  equation  also  has  its  basis  in  the  Boltzmann  equation  —  the  former 
can  be  derived  from  the  latter  through  the  Chapman- Enskog  analysis  (14).  That  very  same 
Chapman-Enskog  analysis  can  be  used  to  show  that  the  lattice  Boltzmann  methodology  can 
be  applied  to  solve  any  conservation  law  of  the  continuous  Boltzmann  equation  including  the 
Navier-Stokes  equations.  It  has  also  been  proved  that  the  lattice  Boltzmann  equation  tan¬ 
tamount  to  an  explicit  finite  difference  scheme  of  the  Navier-Stokes  equations  with  second 
order  spatial  accuracy  and  first  order  temporal  accuracy  with  a  non-zero  compressibility 
(11;  12;  13).  The  present  day  lattice  Boltzmann  equation,  with  its  high-fidelity  physics 
and  computation-efficient  formulation,  is  a  viable  alternative  to  the  continuum  methods  for 
simulating  fluid  flows.  In  fact,  it  can  be  argued  that  for  many  complex  problems  involving 
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multi-fluid  phenomena,  the  physics  can  be  more  naturally  captured  by  the  Boltzmann- 
equation  based  methods  rather  than  Navier-Stokes  equation  based  methods.  Recently  LBE 
method  has  been  extended  to  multiphase  flows  (15;  16;  17;  18)  and  multicomponent  flows 
(19;  20;  21;  22;  23;  24;  25;  26;  27;  28;  29;  30;  31),  flows  through  porous  media  (32;  33), 
and  particulate  suspensions  in  fluids  (34;  35;  36;  37).  Most  existing  LBE  models  for  multi- 
component  fluids  (19;  20;  21;  22;  23;  24;  25;  26;  27)  tend  to  be  somewhat  heuristic  and  make 
the  single-fluid  assumption.  The  single-relaxation-time  or  Bhatnagar-Gross-Krook  (BGK) 
approximation  (38)  is  used  in  most  existing  models  (21;  22;  23;  24;  25;  26;  27)  restricting 
applicability  to  unity  Prandtl  and  Schmidt  numbers. 

A  rigorous  mathematical  development  of  multi-fluid  lattice  Boltzmann  equation  for  mul¬ 
ticomponent  fluids  is  still  in  its  infancy  and  such  is  the  object  of  the  present  work.  As  a  first 
step,  in  this  work  we  develop  a  two-fluid  lattice  Boltzmann  model  which  is  based  on  kinetic 
theory  for  binary  mixtures.  Such  a  model  would  be  capable  of  (i)  simulating  arbitrary 
Schmidt  and  Prandtl  numbers,  and  (ii)  accurately  modeling  the  interaction  between  mis¬ 
cible  and  immiscible  fluids.  We  follow  a  general  approach  within  the  framework  of  kinetic 
theory  for  developing  the  lattice  Boltzmann  models  for  multi-fluid  mixtures.  This  work 
is  a  part  of  our  ongoing  effort  to  set  the  lattice  Boltzmann  equation  on  a  more  rigorous 
theoretical  foundation  and  extend  its  use  to  more  complex  flows.  We  derive  a  discretized 
version  of  the  continuum  Boltzmann  equations  for  binary  mixtures.  The  extension  of  this 
methodology  to  multi-fluid  mixtures  is  relatively  straight-forward. 

Kinetic  theory  of  gas  mixtures  has  received  much  attention  in  literature  (14;  39;  40;  41; 
42;  43;  44;  45;  47;  46;  48;  49;  50;  51).  Many  of  the  kinetic  models  for  gas  mixtures  are 
based  on  the  linearized  Boltzmann  equation  (38;  52;  53),  especially  the  single-relaxation¬ 
time  model  due  to  Bhatnagar,  Gross,  and  Krook  —  the  celebrated  BGK  model  (38).  The 

kinetic-theory  mixtures  model  employed  in  this  work  was  proposed  by  Sirrovich  (43)  which 
is  also  linear  in  nature. 

This  document  is  a  detailed  follow-up  to  our  previous  work  published  as  a  Rapid  Com- 
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munications  in  PRE  (31)  and  it  is  organized  as  follows.  Section  B.2  provides  a  brief  review 
for  some  of  the  existing  kinetic  models  of  mixtures  that  form  the  theoretical  basis  of  the 
present  work.  Section  B.3  contains  the  derivation  of  the  lattice  Boltzmann  model  for  binary 
mixtures  from  the  corresponding  continuous  Boltzmann  equations.  In  Sec.  B.4,  the  hydro- 
dynamic  equations  of  the  lattice  Boltzmann  model  are  determined.  Section  B.5  contains  the 
derivation  of  the  diffusion  force  and  the  mutual  diffusion  coefficient  in  the  lattice  Boltzmann 
model,  and  the  diffusion-advection  equations  for  the  mass  and  molar  concentrations  of  the 
system.  Section  B.6  discusses  the  short  and  long  time  behaviors  of  the  model.  Section  B.7 
concludes  the  Appendix  with  a  summary  of  the  present  work  and  possible  directions  of 
future  work.  The  three  appendixes  contain  the  details  of:  B.A.l  the  iterative  procedure  to 
solve  the  Boltzmann  equation;  B.A.2  the  discretized  equilibrium  distribution  function;  and, 
B.A.3  the  Chapman-Enskog  analysis  of  the  lattice  Boltzmann  model  for  binary  mixtures. 

B.2  Kinetic  theory  of  gas  mixtures 

Following  a  procedure  similar  to  the  derivation  of  the  Boltzmann  equation  for  a  pure  system 
of  single  species,  one  can  derive  N  simultaneous  equations  for  a  system  of  N  species  by 
reducing  the  appropriate  Liouville  equation.  For  the  sake  of  simplicity  without  loss  of 
generality,  we  shall  only  discuss  the  Boltzmann  equations  for  a  binary  system  here.  The 
simultaneous  Boltzmann  equations  for  a  binary  system  are 

dtfA  +  £■' V/A  +  aA- V?/A  =  Qaa  +  Qab,  (B.la) 

dtfB  +  £-V/B  +  oB -Ve/B  =  Qba  +  Qbb,  (B.lb) 

where  QAB  =  QBA  is  the  collision  term  due  to  the  interaction  among  two  different  species 

A  and  B.  Obviously,  for  an  A"- component  system,  there  will  be  N  such  equations,  each 

containing  N  collision  terms  on  the  right  hand  side.  In  general,  the  collision  term  is 


Qab  =  I d&<m<TABuB  -  all  [fAf'B  -  fAfB] 


(B.2) 
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where  ctab  is  the  differential  cross  section  of  A-B  collision,  and  Q  is  the  solid  angle,  and 
fiA  (y'B)  an(j  jA  (yB)  denote  the  post-collision  and  pre-collision  state  of  the  particle  A 
(B),  respectively.  Obviously,  the  equations  for  a  system  of  multiple  species  are  much  more 
formidable  to  analyze  than  the  comparable  Boltzmann  equation  for  a  pure  system  of  single 
species.  The  first  modeling  objective  is  to  find  a  suitable  approximation  for  the  collision  term 
given  by  Eq.  (B.2)  that  would  substantially  simplify  computation  without  compromising 
the  essential  physics. 

The  justified  approximation  of  the  collision  terms  must  rely  on  our  clear  understanding 
of  the  underlying  physical  system.  In  a  system  of  multiple  species,  there  are  a  number  of 
competing  equilibration  processes  occurring  simultaneously.  The  approach  to  equilibrium 
in  the  system  can  be  roughly  divided  into  two  stages.  First  of  all,  each  individual  species 
equilibrates  within  itself  so  that  its  local  distribution  function  approaches  a  local  Maxwellian 
distribution.  This  process  of  individual  equilibration  is  also  referred  to  as  Maxwellization. 
In  the  second  stage,  the  entire  system  equilibrates  so  that  the  velocity  and  temperature  dif¬ 
ferences  among  different  species  vanishes  eventually.  There  are  many  different  time  scales  in 
the  equilibrating  process  of  a  multicomponent  system.  In  addition,  the  Maxwellization  itself 
can  take  place  in  many  scenarios  depending  on  the  molecular  weights  and  mass  fractions  of 
the  participating  species.  Consider  two  mixtures,  each  consisting  of  a  light  and  heavy  gas. 
In  Mixture  1,  the  total  mass  of  each  species  is  the  same,  implying  smaller  number  density 
for  the  heavier  gas.  In  Mixture  2,  the  number  densities  of  the  two  species  is  the  same, 
implying  the  mass  density  (or  mass  fraction)  of  the  heavier  species  is  larger.  In  Mixture  1, 
the  Maxwellization  of  light  species  is  mostly  due  to  self-collision  whereas  the  equilibration  of 
the  heavier  species  is  predominantly  due  to  cross  collisions.  This  is  due  to  the  fact  that  the 
number  of  molecules  of  heavy  species  available  for  collisions  is  small.  In  Mixture  2,  where 
the  number  densities  of  the  two  species  are  comparable,  Maxwellization  of  both  species  in¬ 
volves  self  and  cross  collisions.  When  the  process  of  Maxwellization  is  complete,  the  stress 
of  the  corresponding  species  becomes  isotropic,  or  equivalently  the  heat  conduction  relaxes. 
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Therefore,  the  scale  on  which  the  stress  becomes  isotropic  or  the  heat  conduction  relaxes  is 
a  suitable  measure  of  Maxwellization. 

The  equilibration  among  different  species  can  also  take  place  in  several  different  manners. 
Velocity  and  temperature  differences  may  equilibrate  on  the  same  temporal  scale  (as  in  Mix¬ 
ture  1  above)  or  on  vastly  different  scales  (as  in  Mixture  2).  In  addition,  these  equilibrating 
processes  need  not  to  occur  sequentially  but  also  concurrently  with  the  Maxwellization. 

There  is  a  significant  amount  of  literature  on  gas  mixtures  within  the  framework  of 
kinetic  theory  (14;  39;  40;  41;  42;  43;  44;  45;  47;  46;  48;  49;  50;  51).  In  the  Chapman- 
Enskog  analysis  for  a  simple  gas,  one  assumes  a  clear  separation  of  scales  in  space  and 
time,  that  is,  to  distinguish  the  spatial  and  temporal  scales  which  are  much  larger  than 
the  mean  free  path  or  mean  free  time.  An  analogy  for  a  mixture  becomes  difficult  because 
of  multiplicity  of  length  scales.  In  the  classic  work  of  Chapman  and  Cowling  (14),  the 
full  Boltzmann  equations  (with  integral  collision  terms)  for  a  binary  mixture  are  analyzed 
under  the  assumptions  that  all  scales  are  roughly  of  the  same  order,  or  equivalently,  that  the 
phenomenon  to  be  examined  is  smooth  with  respect  to  all  collisional  scales.  Determination 
of  the  various  transport  coefficients  —  viscosities,  diffusivities,  thermal  diffusivities  and 
conductivity  —  was  the  main  objective  of  that  work.  However,  no  attempt  was  made  to 
describe  the  dynamics  of  the  evolution. 

Direct  analysis  or  computation  of  the  Boltzmann  equation  is  not  generally  feasible. 
This  is  due  to  the  difficulty  involved  in  evaluating  the  complex  integral  collision  operators. 
To  make  further  progress  one  can  follow  one  of  two  approaches.  The  first,  Grad’s  moment 
method,  is  to  obtain  the  non-normal  solutions  of  the  Boltzmann  equation  (i.e.,  the  solutions 
beyond  the  hydrodynamic  or  conserved  variables)  (54).  Closure  modeling  would  then  be 
required  to  express  the  unclosed  moments  in  terms  of  the  closed  moments.  And  the  second 
is  to  derive  simplified  model  equations  from  the  Boltzmann  equation  which  are  more  man¬ 
ageable  to  solve.  Many  model  equations  are  influenced  by  Maxwell’s  approach  to  solving  the 
Boltzmann  equation  by  making  extensive  use  of  the  properties  of  the  Maxwell  molecule  (55) 
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and  the  linearized  Boltzmann  equation.  The  simplest  model  equations  for  a  binary  mixture 
is  that  by  Gross  and  Krook  (41),  which  is  an  extension  of  the  single-relaxation- time  model 
for  a  pure  system  —  the  celebrated  BGK  model  (38). 

With  the  BGK  approximation  (38;  41),  the  collision  integrals  [a,  q  6  (A,  B)]  can 
be  approximated  by  following  linearized  collision  terms 

(B.3a) 

=  ~ir  -  <B.3b) 

where  and  /<T?(°)  are  Maxwellians 


M°)  =  n<r  p-(£~u<t)2/(2R<tT',) 

(. 2irRaTa)D /2 

°)  =  _ n<r  r—(£—u„c)2/(2  R„T„A 

(27 xR^)0'2 


(B.4a) 

(B.4b) 


where  D  is  the  spatial  dimension,  Ra  =  kB/ma  is  the  gas  constant  of  the  a  species,  kB  is 
the  Boltzmann  constant  and  ma  is  the  molecular  mass  of  the  a  species.  There  are  three 
adjustable  relaxation  parameters  in  the  collision  terms:  A^,  Af,  and  Aw  =  {njna)\a. 
The  first  Maxwellian  is  characterized  by  the  conserved  variables  of  each  individual 
species:  the  number  density  nCT,  the  mass  velocity  ua,  and  the  temperature  Ta;  while 
the  second  Maxwellian  /^(°)  and  /?<T(°)  is  characterized  by  four  adjustable  parameters: 
u<rqt  u^cr,  Taq,  and  T^fj.  There  are  several  considerations  in  determining  these  arbitrary 
parameters:  simplicity  of  the  resulting  theory,  accuracy  of  approximation,  and  ease  of 
computation.  Cross-collisional  terms  will  be  symmetric  only  if  one  takes  uaq  =  it?cr  =  u 
and  Ta<;  =  TCiCr  =  T ,  where  u  and  T  are  the  velocity  and  temperature  of  the  mixture.  This  is 
essential  in  preserving  a  similarity  to  irreversible  thermodynamics,  especially  the  Onsager 
relation  (56).  On  the  other  hand,  fewer  terms  in  the  expansion  of  fa  about  /CTC(°)  would 
be  needed  in  many  cases  if  one  chooses  ua<;  =  ua  and  =  Tff,  i.e.,  /a?(°)  =  /CT(°).  One 
salient  difference  between  using  u  and  T  of  the  mixture  in  the  Maxwellian  as  opposed 
to  using  ua  and  Ta  for  the  species  is  that  the  former  choice  leads  to  a  single-fluid  theory, 
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i.e.,  a  set  of  hydrodynamic  equations  for  the  mixture,  while  the  latter  leads  to  a  two-fluid 
theory  (43;  49),  i.e.,  two  sets  of  hydrodynamic  equations  for  the  species.  Obviously,  in  the 
cases  where  the  properties  between  the  two  species  are  vastly  different,  the  two-fluid  theory 
is  preferred  (56). 

The  cross-collision  term  Jaq  can  be  better  approximated  by  expanding  /<T  around  the 
Maxwellian  (43) 


f<*(  o) 

Tlo-h  bTq- 


Pd^c  '  (Wff  ^?) 


3 

+/iT2 

~Ma, 


(-£- 

\2RcT,, 

(— 

\2RaTa 
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where  ca  =  (£  —  ua)  is  the  peculiar  (or  thermal)  velocity  of  the  a  species,  and 


(B.5) 


—  pm 


P&Ps 

P 


_ |_  „/  (n<r  7h) 

pa  mn(Tn?(m0.  +  m?)J  ’ 


(B.6) 


and  pD,  pt .  pm,  p'm  are  positive  and  at  most  functions  of  density  and  temperature  (43), 
the  physical  significance  of  these  parameters  are  to  be  discussed  next. 

We  now  consider  the  following  model  equations  for  a  binary  mixture  due  to  Sirovich 
(43): 


dtfff  +  Z-Vr  +  a<,'Vtfr  =  J**  +  J",  (B.7) 


where  the  self-collision  term  Jaa  is  approximated  with  the  BGK  model  of  Eq.  (B.3a),  and 
the  cross-collision  term  J(T?  is  given  by  Eq.  (B.5).  When  the  external  force  is  not  present 
(acr  =  0),  and  if  p,D ,  //-t ,  Pm ,  p'm  are  considered  as  constants  (for  the  sake  of  simplicity)  we 
can  immediately  obtain  the  following  moment  equations  from  the  above  equations 


dt(Ua  W^)  —  Pd 


(ua  ti?), 


dt(Ta-T<)  =  -pT±- 

+JL(Mzi 

3 kB  V  n<r 


(—  +  — )  (Ta  ~  T<) 
\ria  njx 

nf  )  {U<r  ~  U ^  ■ 


(B.8a) 


(B.8b) 
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The  above  equations  describe  the  exponential  decay  of  the  velocity  and  temperature  dif¬ 
ferences  for  the  two  species,  as  discussed  in  the  earlier  remarks  regarding  the  processes  of 
Maxwellianization  and  equilibration  in  the  mixture.  The  physical  significance  of  the  pa¬ 
rameters  /j,D,  pT,  Mm,  and  n'm  become  apparent  in  the  above  equations  —  these  parameters 
determine  the  relaxation  rates  in  the  Maxwellization  processes. 

Solving  Eqs.  (B.7)  by  means  of  iteration  (cf.  Ref.  (43)  or  Appendix  B.A.l  ,  one  first 
obtains 


Ucr  —  —  U, 

Ta  —T^  =  T, 

r0)  =  /«(o)(fVr>  T) 


n0 


«-«) 


21 


(B.9a) 

(B.9b) 

(B.9c) 


{2ttRi7T)d/2  6XP  L  2 R„T 

Substituting  the  above  results  into  the  left-hand  side  of  Eqs.  (B.7),  one  has  the  following 
equation  for  the  second-order  solution  of  fa  (43): 


rt(0)(»„  u,T)\-C'-d„  + 

L  n<r 

+(-£--?) 

\2RaT  2) 


C(ji  C(jj 

~Rjf 
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.  —  [f<T  -  H0)]  - 
KU  J  1  naksTc 
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where  cc  =  (£  -  ua),  and  the  diffusion  force  da  and  the  rate-of-shear  tensor  are  given 
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by 
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ij'j  1 


(B.lla) 

(B.llb) 


and  p  =  nkBT  is  the  total  pressure  of  the  mixture.  Prom  the  above  solution  of  /CT,  one 
can  compute  for  the  relative  velocity  {ua  -  uj,  the  temperature  difference  (Ta  -  Tf),  the 
(traceless)  stress  tensor  pij,  and  the  heat  flux  q  (43): 


nkBT  ,  n“ 


<  =  - 


{uc  u^)  — 

Pd 

_T1_  PoKn  -Tp 


(Ta-Tq)  = 


Da,d'a, 
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q  =  -&bT [ncr(u<7  -  u)  +  «?(«?  -  «)] 

-hT{t+\)kBVT’ 


where  the  diffusion  force 


PaP<, 

pp 


T  np 

(ciu  —  a?), 


and  the  binary  diffusion  coefficient  (53)  of  the  model 


jj - - 

npD 


(B.12a) 

(B.12b) 

(B.12c) 

(B.12d) 


(B.13a) 


(B.14) 


The  self-diffusion  coefficient  D aa  is  a  special  case  of  the  above  formula  when  n  =  n?  =  na. 
The  viscosity  v  and  thermal  conductivity  k  can  also  be  read  from  the  above  formulas  for 
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Pij  and  q  as  the  following 

v  —  ksT  >  (B.15a) 

-Hr(£  +  £)-  <B-15b> 

The  above  transport  coefficients  are  determined  by  the  parameters  A^  and  pD,  //r,  pm, 
and  p'm:  Xa  determines  the  viscosity  and  the  thermal  conductivity  of  the  a  species  and  the 
combination  of  Ao-’s  determines  that  of  the  mixture;  pD  determines  the  diffusion  coefficients 
in  the  model;  and  fJ.Dp'n ljpT  determines  the  diffusion  of  the  temperature  difference  due  to 
velocity  difference. 

Two  salient  features  of  the  model  described  by  Eqs.  (B.7)  should  be  addressed.  First, 
the  cross-collision  term  J<T?  of  Eq.  (B.5)  is  exact  for  the  Maxwell  molecules  obeying  the 
inverse  fifth-power  interaction  potential.  Equations  (B.7),  therefore,  can  be  considered  to 
be  a  model  for  the  Maxwell  gas  (43).  One  immediate  consequence  of  this  approximation  is 
that  the  diffusion  force  of  Eq.  (B.12a)  does  not  contain  a  thermal  diffusion  term,  as  it  should. 
Second,  the  BGK  approximation  of  the  self-collision  term  Ja(I  of  Eq.  (B.3a)  imposes  the 
limitation  of  a  unity  Prandtl  number.  However  both  these  limitations  of  the  model  can  be 
overcome  by  using  the  linearized  Boltzmann  equation  (57;  58;  59)  with  multiple  relaxation 
times  and  a  nonlinear  approximation  of  the  collision  terms  (43;  49;  60). 

B.3  The  lattice  Boltzmann  model  for  binary  mixture 

We  shall  construct  a  lattice  Boltzmann  model  for  binary  mixtures  based  on  the  model 
given  by  Eqs.  (B.7).  In  the  present  work  we  only  consider  the  isothermal  case  such  that 
Tar  =  Tf  =  Taq  =  T  =  const.  Consequently,  we  can  also  ignore  the  terms  related  to  thermal 
effects  in  Ja<;  of  Eq.  (B.5)  by  setting  pT  =  Mm  =  p'm  =  0,  i.e., 


1  Or 


(B.16) 
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where  the  equilibrium  function  for  the  cr  species  is  chosen  to  be  the  Maxwellian 

equilibrium  distribution  depending  on  the  mass  velocity  of  the  a  species  ua  as  the  following 


0)  = 


Per 


exp 
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2 RaT  J 
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Note  that  from  hereafter  and  fa  are  the  single  particle  mass  density  distribution 
functions,  as  opposed  to  the  single  particle  number  density  distribution  functions.  We  can 
derive  the  lattice  Boltzmann  equation  by  discretizing  the  model  equations  (B.7),  following 
the  procedure  described  in  (4;  5): 


fa  (®i  +  eafit,t  +  $t)  ~ 

=  Ja+J«-F°Ju 


(B.18) 


where  the  self-collision  term  J™,  the  cross-collision  term  ,  and  the  forcing  term  Fa  are 
given  by 


K’  =  -  “)•(“»  “  «t). 

td  p  Ral 

t pa  6a-Q,0- 

Fa  =  ~WaP<r 


(B.19a) 

(B.19b) 

(B.19c) 


RaT  ’ 

where  pa  and  and  ua  and  are  the  mass  densities  and  flow  velocities  for  species  cr  and 
<;,  they  are  the  moments  of  the  distribution  functions: 

p»  =  E^  =  E^<0).  <B-20a> 

a  a 

PaUa  =  E  =  E  ^(°)e“’  (B.20b) 

Ot  Q 

and  p  and  u  are  respectively  the  mass  density  and  the  barycentric  velocity  of  the  mixture: 


P  “  Pa  +  Pqi 

p\L  —  P(jXL(j  +  pqUq. 
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The  collision  terms  J£a  and  J°q  are  constructed  in  such  a  way  to  respect  the  mass  and 
momentum  conservation  laws.  (The  derivation  of  the  forcing  term  F£  is  given  in  Refs.  (17; 
18).) 

The  equilibrium  distribution  function  /q has  the  following  form  in  general  (cf.  Ap¬ 
pendix  B.A.2  ): 
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where  coefficients  {iaQ}  depend  on  the  discrete  velocity  set  {eQ}.  For  the  sake  of  concrete¬ 
ness  and  simplicity  without  losing  generality,  we  shall  restrict  ourselves  to  a  nine-velocity 
model  on  a  two-dimensional  square  lattice  (D2Q9  model).  In  this  case, 


4/9, 
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1/9, 
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(B.23) 

1/36, 

a  =  5  -  8. 

B  .4  Hydro  dynamics 


The  left-hand  side  of  Eq.  (B.18)  can  be  expanded  in  a  Taylor  series  in  St  (up  to  second 


order  in  St)  and  the  equation  can  be  rewritten  as: 

StDafZ  +  \s\Dlfl  =  J r  +  J?  -  FZSt,  (B.24) 

where  Da  =  dt  +  ea  •  V.  Obviously, 

E-,r  =  E^'  =  E^=°.  (B-25a> 
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By  means  of  the  Chapman-Enskog  analysis  (multiple-scale  expansion),  we  can  derive  the 
hydrodynamic  equations  for  the  mixture  from  Eq.  (B.24)  (see  details  in  Appendix  B.A.3  ). 

The  mass  conservation  laws  for  each  species  and  the  mixture  can  be  derived  immediately 
from  Eq.  (B.24): 
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Note  that  the  mass  conservation  does  not  hold  for  each  individual  species  at  the  Navier- 
Stokes  level,  although  it  does  at  the  Euler  level.  However,  the  mass  conservation  law  does 
apply  to  the  mixture  as  a  whole.  The  right  hand  side  of  Eq.  (B.26)  reflects  the  mass  flux 
due  to  diffusion. 

We  can  also  derive  the  Euler  equation  for  each  species: 


Pad^ll,?  -f-  PpUtf'^lla  —  V pa 

+p<TaCT  -  -j- Pa*  (ua  -  ug),  (B.28) 

where  pa  =  naksT  =  pcRaT  is  the  partial  pressure  of  the  a  species,  and  the  Navier-Stokes 
equation: 
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where  the  viscosity  of  the  cr  species  is 
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Eq.  (B.29)  is  consistent  with  the  results  in  Ref.  (49). 
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B.5  Diffusion  in  isothermal  mixtures 


The  difference  between  the  two  Navier-Stokes  equations  for  individual  species  (cr  and  <r) 
leads  to  the  following  equation: 


-y  K  -  ««)  =  —~rd^ 

PcrPs 

—  {dt&u  +  'V-(uSu)  +  V2(i —  i/5us)}  , 


(B.31) 


where  6u  =  (tv  —  uq ),  u  =  \{ua  +  tic),  X?-(u5u)  =  u-'VSu  +  <5tx- Vtt,  and  the  diffusion 
force 


=  v(lf)  +  (^-7)vtap 

+^(a?  -  a„)  (B.32) 

PP 

=  V  (?)  VlnP 

\  n  /  np 

+  __  (aC  ~  acr) 

PP 

~  — 


where  p  =  nfc^T  is  the  total  pressure  of  the  mixture,  and  the  total  number  density  n  is 

n  =  nCT +nf  =  —  +  — .  (B.33) 

ma  m? 

The  diffusion  force  includes  the  effects  due  to  the  molar  concentration  gradient  V(n(7/n),  the 
total  pressure  gradient  and  the  particle  mass  difference  (ma  —  mf) V  lnp,  and  the  external 
force  ( aa  —  aq). 

It  has  already  been  assumed  in  the  derivation  of  the  two-fluid  equations  that  derivatives 
are  slowly  varying  on  the  time  scale  of  Maxwellization  (49).  Therefore  the  terms  inside  the 
curly  brackets  Eq.  (B.31)  can  be  neglected  in  the  diffusion  time  scale.  Thus,  to  the  leading 
order,  we  have 

(tV  -  tv)  =  -TDbt-^—da.  (B.34) 

PcrPs 
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Also,  by  definition  (53,  Eq.  (6.5-7a)),  we  have 


n 


(tier  t l,)  —  Dffqdff, 

TtcfTlq 

thus  the  mutual  diffusion  coefficient  in  the  mixture  is 


(B.35) 


(B.36) 


Da,  =  ~~~~~TD8f 
nrriam c 

Note  that  the  difference  between  the  above  formula  and  Eq.  (B.14)  is  due  to  the  difference 
of  a  factor  Pap,/p  between  of  Eq.  (B.19b)  and  Jaq  of  Eq.  (B.5),  i.e., 

P  Pd 

The  mass  flux  of  the  a  species,  by  definition,  is 


3*  =  Pa{Ua  ~  U)  =  -TD5tpda, 


where  we  have  used  the  identity  that  p(ua  —  u)  =  p,(ua  —  uc).  The  continuity  equa¬ 
tion  (B.26)  for  the  a  species  can  be  written  as 


Dtpa  +  P<rV-U  +  V-ja  =  —  2^*^  (pda)  ,  (B.37) 

where  Dt  =  (dt  +  u-V).  By  assuming  the  incompressibility  of  the  fluid  (i.e.,  V  ti  =  0),  we 
obtain  the  following  advection-diffusion  equation  for  an  isothermal  mixture 

dtpa  +  u-Vpa  =  V-T*Stpda,  t*  =  (td  -  0  .  (B.38) 

Similar  to  the  Henon  correction  for  the  viscosity  (61),  the  diffusivity  is  modified  by  the 
second-order  discrete  effect: 


_  PksT  (  _  1\  .  _  pksT 

nmaTn q  \  D  2 )  1  nmam  D  1 


(B.39) 


Obviously,  the  self-diffusion  coefficient  in  the  lattice  Boltzmann  model  for  the  a  species  is 
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The  mass  concentration  <fi  and  molar  concentration  <p  (dimensionless  order  parameters) 


are  defined  as 

,  [pc  Z  EA  m  -  ~  n?) 

{pc  +  Pc)  ’  (tv  +  nc)  ’ 

and  they  related  to  each  other 

_  (m^  -  mc)  +  {jug  +  rn^yp 
{' ma  +  mc)  4-  {ma  -  m?)<p  ’ 

{mq  -  ma)  +  (mc  +  ma)(f) 

(p  =  - # 

(mq  +  ma)  +  (mq  -  ma)<f>' 


(B.41) 


(B.42a) 

(B.42b) 


The  diffusion  force  can  be  written  in  terms  of  (f>  and  <p: 


V<p  +  (<p  -  4>)W  lnp  +  ^  ~ (Qc  -  «<r) 


The  nonlinear  diffusion-advection  equations  satisfied  by  </>  and  </?  can  be  derived  from 
Eq.  (B.38): 

+  u  • ' V<t>  =  i' V  •  {D%V4>  +  F) ,  (B.43a) 

dtip  +  u-'V<p  =  ^V-{D*vV<p  +  F),  (B.43b) 

where 


pr’DS,  //V  kBT  ,/‘ 


D.  = 

v  mam^  \n/ 

D;=pr*St 

F  = 


n 


r*DSt, 


(<P  ~  4>)^P  +  2 Pi1  ~  02)(°c  -  °c ) 
l1  ~<P)  +  (1  +  <P) 


r*DSt, 


m„ 


mc 


(B.44a) 

(B.44b) 

(B.44c) 

(B.44d) 


Obviously,  when  ma  =  mc,  <p  =  <f>,  and  then  Eq.  (B.43a)  and  Eq.  (B.43b)  become  identical. 


B.6  Short  and  long  time  behaviors 

As  discussed  in  Sec.  B.2,  short  and  long  time  behaviors  of  a  binary  mixture  involve  different 
Maxwellization  and  equilibration  processes.  This  is  reflected  in  the  macroscopic  equations 
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derived  in  preceding  sections.  In  an  initial  stage  (short  time),  the  diffusion  velocity  (n(7-uc) 
is  significant,  thus  the  system  is  described  2 (D  +  1)  equations,  i.e.,  two  sets  of  mass  and 
momentum  conservation  laws  for  each  species  given  by  Eqs.  (B.26)  and  (B.29).  As  the 
system  equilibrates  so  that  the  diffusion  velocity  ( u a  —  is  diminishing,  the  physics  is 

then  described  by  (D  +  2)  equations:  the  continuity  equation  of  the  mixture,  Eq.  (B.27), 
the  Navier-Stokes  equation  for  the  barycentric  velocity  of  the  mixture  u,  and  the  diffusion- 
advection  equation  for  the  mass  [Eq.  (B.43a)]  or  molar  [Eq.  (B.43b)]  concentrations.  Only 
in  very  late  stage  of  equilibration,  the  concentration  behaves  more  or  less  as  a  passive  scalar. 

The  Navier-Stokes  equation  for  u  is 

pdtu  +  pu-Vu  =  -Vp  +  V-n  +  pa,  (B.45) 

where  p  =  k,BT(na  +  nc),  pa  =  P(jaa  +  pqaq,  and 

n  =  E./W(VM  +  (v«ff)t] 

~  {p<tV<j  +ft^)[(V?i)  +  (Vu)*].  (B.46) 

In  the  derivation  of  Eq.  (B.45),  we  ignore  two  terms:  one  is  (rDD*?n/i20-Tn0-nc)V(pd0-)2, 
and  the  other  is  in  proportion  of  u(uff  -  ttc)  due  to  Jg*  [cf.  Eq.  (B.58)].  These  terms  are 
negligible  when  the  mixture  is  more  or  less  homogeneous.  Also,  when  the  diffusion  velocity 
(ua  —  uq)  vanishes,  the  approximation  in  Eq.  (B.46)  becomes  exact. 

Some  remarks  to  place  the  present  work  in  perspective  are  in  order.  Unlike  most  existing 
lattice  Boltzmann  models  for  binary  mixtures  (22;  23;  24;  25;  26;  27),  the  mutual  and  self 
diffusion  coefficients  of  the  present  model  are  independent  of  the  viscosity.  (We  note  that  an 
existing  model  proposed  by  Flekkpy  (20)  already  has  this  feature.)  The  diffusion  coeflicients 
depend  on  the  relaxation  parameter  td  and  relevant  physical  properties  of  the  mixture,  such 
as  the  molecular  masses  of  each  species  ma  and  mq,  etc.  The  present  model  is  therefore 
capable  of  incorporating  more  general  physics.  In  addition,  the  present  model  can  simulate 
both  miscible  and  immiscible  binary  mixtures  by  changing  the  sign  of  (rD  -  1/2),  i.e.,  for 
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positive  ( To  -  1/2),  the  mixture  is  miscible,  and  for  negative  (rD  -  1/2),  the  mixture  is 
immiscible.  We  have  performed  numerical  simulations  to  verify  that:  (a)  the  diffusion 
coefficient  D*q  depends  on  td  as  given  by  Eq.  (B.39)  and  it  is  independent  of  ra  and  t?; 
and  (b)  when  ( td  — 1/2)  <  0,  spinodal  decomposition  or  separation  (anti-diffusion)  between 
different  species  in  the  mixture  occurs.  (The  details  of  the  numerical  results  will  be  reported 
elsewhere.) 

B.T  Discussion  and  conclusions 

We  have  constructed  a  lattice  Boltzmann  model  for  binary  mixtures  with  several  important 
features.  All  the  modeling  issues  are  addressed  at  the  continuum  level  within  the  framework 
of  extended  kinetic  theory.  The  lattice  Boltzmann  model  is  then  directly  derived  from  the 
continuous  kinetic  model  equations  using  a  formal  discretization  procedure.  The  lattice 
model  thus  inherits  the  sound  physics  and  mathematical  rigor  incumbent  in  kinetic  theory. 
This  is  in  contrast  to  previous  lattice  Boltzmann  models  for  mixtures  (22;  23;  24;  25;  26;  27), 
which  are  not  directly  based  on  the  fundamental  physics  of  continuum  kinetic  equations. 
These  models  rely  on  fictitious  interactions  (22;  23)  or  heuristic  free  energies  (24;  25;  26;  27) 
to  produce  the  requisite  mixing.  (Many  defects  of  the  free-energy  models  (24;  25;  26;  27)  are 
due  to  the  incorrectly  defined  equilibria  (18).)  These  non-physical  effects  present  a  further 
problem  since  they  are  not  easily  amenable  to  mathematical  analysis  (17;  18).  The  heuris¬ 
tic  elements  of  the  previous  lattice  Boltzmann  models  (22;  23;  24;  25;  26;  27)  have  been 
eliminated,  resulting  in  a  physically  justifiable  model  that  is  simple  to  compute.  Further, 
due  to  the  close  connection  to  kinetic  theory,  the  derivation  of  the  hydrodynamic  equations 
associated  with  the  lattice  Boltzmann  model  is  significantly  simplified  and  rendered  math¬ 
ematically  more  rigorous.  The  derivation  of  the  hydrodynamic  equations  from  the  previous 
lattice  models  are  much  less  rigorous  (22;  23;  24;  25;  26;  27). 

The  second  important  feature  of  the  present  work  is  that  the  model  is  based  upon  a  two- 
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fluid  theory  of  binary  mixtures.  The  previous  models  (19;  20;  21;  22;  23;  24;  25;  26;  27;  51), 
on  the  other  hand,  are  derived  from  a  simpler,  but  highly  restrictive,  one-fluid  theory. 
In  the  single-fluid  models  with  BGK  approximation  one  is  constrained  to  use  the  ad  hoc 
“equilibrium  velocity”  (22;  23;  51) 

u(eq )  _  ('frPtr^o-  +  Tg^U^) 

(T?/V  +  TcrPc) 

in  the  equilibrium  distribution  function  /q(0)  in  order  to  satisfy  the  local  conservation  laws. 
As  a  result,  the  viscous  relaxation  process  and  the  diffusion  process  are  inseparable.  The 
analysis  of  these  models  therefore  becomes  unnecessarily  tedious  and  cumbersome  (22;  23). 
The  models  with  free  energies  (24;  25;  26;  27)  do  not  yield  correct  hydrodynamic  equations 
(17;  18),  mostly  due  to  the  incorrectly  defined  equilibrium  distribution  functions  used  in 
these  models  (18;  62).  Furthermore,  single-fluid  models  cannot  be  applied  to  mixtures  of 
species  with  vastly  different  properties.  In  the  present  two-fluid  model,  the  diffusion  behav¬ 
ior  is  decoupled  from  viscous  relaxation.  The  diffusivity  is  determined  by  the  parameter 
td  and  the  physical  properties  of  the  mixture.  The  model  is  capable  of  simulating  either 
miscible  or  immiscible  fluids  by  changing  the  sign  of  (td  —  1/2). 

The  proposed  LBE  model  for  binary  mixtures  simulates  diffusion  by  considering  a  mu¬ 
tual  interaction  term  leading  to  the  diffusion  velocity  (ua  -  u?)  which  is  directly  related  to 
the  diffusion  driving  force  in  binary  mixtures.  The  diffusion  velocity  (up  —  wc)  is  of  the  first 
order  in  the  density  gradient  V p.  This  suggests  that  the  proposed  model,  however,  does 
not  include  any  higher  order  terms  of  the  density  gradient.  This  in  turn  implies  that  the 
proposed  model  does  not  have  a  surface  tension,  which  is  related  to  the  density  gradient 
square  |  Vp|2.  To  include  the  effect  due  surface  tension,  the  terms  related  to  |Vp|2  must  be 
explicitly  considered. 

To  further  improve  the  proposed  lattice  Boltzmann  model  for  binary  mixing,  efforts 
are  currently  underway  are:  (a)  development  of  a  multiple-relaxation-time  model  for  the 
self-collision  term  (57;  58;  59)  that  will  significantly  enhance  the  numerical  stability  of  the 
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scheme  (58;  59);  (b)  consideration  of  models  with  surface  tension;  (c)  inclusion  of  thermal 
diffusion  effects  which  may  be  important  in  combustion  applications;  and,  (d)  development 
of  a  model  for  non-Maxwellian  molecules.  It  should  be  noted  that  all  the  existing  lattice 
Boltzmann  models  are  only  applicable  to  the  Maxwell  molecules  as  a  direct  consequence 
of  the  linearization  of  the  Boltzmann  equation.  This  limitation  can  be  overcome  by  either 
using  a  different  expansion  of  the  distribution  function,  or  by  including  the  non-Maxwellian 
effects  in  the  collision  terms. 


Appendix 

B.A.l  Iterative  Solution  of  the  Boltzmann  Equation 


We  present  a  short  description  of  the  iterative  procedure  (43)  to  solve  the  Boltzmann 
equation 

Dtf  =  Q[f]  (B.47) 

where  Dt  =  dt  +  £-V  and  Q[f]  is  the  collision  term.  The  (n  +  l)-th  iteration  solution  /(n+1) 
is  obtained  from  the  n-th  iteration  solution  / by  solving  the  equation 


Q[f(n+ 1)]  =  A/(n), 


(B.48) 


subject  to  the  following  conservation  constraints 


(B.49) 


where  /W  =  /(°)(p(7l+1)J  u(n+1\  T(n+1))  is  the  equilibrium  distribution  function  which  de¬ 
pends  on  the  hydrodynamical  moments  p,  n,  and  T  computed  from  y(n+1).  Following  the 
Chapman-Enskog  procedure,  the  temporal  derivatives  are  removed  by  using  the  conserva- 
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tion  equations: 

dtP{n+1)  =  -V-(pu), 

dtu(n+1)  =  -ti.Vti  +  o  -  -  V-P(n), 

P 

dte(n+1)  =  -V-(eu)  +  pa-u  -  V-g(n)  +  P<B> :  Vu, 

where  e  is  the  internal  energy  and  e  =  ynfcsT  for  ideal  gases,  and  P  and  q  are  the  stress 
tenor  and  the  heat  flux,  respectively.  On  right  hand  side  of  the  above  equations,  the 
hydrodynamical  moments  p,  u,  and  e  are  computed  from  the  (n  +  l)-th  iteration  solution 
/(n+1)  but  their  superscript  (n  + 1)  are  omitted  since  they  are  conserved  quantities.  On  the 
other  hand,  the  stress  tensor  P  and  the  heat  flux  g,  which  axe  not  conserved,  are  denoted 
with  the  superscript  n  as  they  are  obtained  from  the  solution  of  the  previous  iteration  f(n\ 
In  general,  the  iterative  procedure  described  above  is  expected  to  converge  more  rapidly 
than  a  procedure  of  successive  approximation,  such  as  the  Chapman-Enskog  procedure 
for  the  Boltzmann  equation.  The  reason  is  that  in  the  iterative  procedure  the  (nonlin¬ 
ear)  integral  equation  must  be  solved  at  each  step  as  given  by  Eq.  (B.48),  whereas  in  the 
Chapman-Enskog  procedure  the  integral  equation  is  only  solved  at  the  initial  step  and  at 
all  approximation  of  higher  order  only  the  linearized  integral  equation  is  solved. 


B.A.2  The  equilibrium  distribution  function 


We  consider  the  equilibrium  distribution  function  of  Eq.  (B.17)  for  a  species  based  on  its 
mass  velocity  ua.  The  distribution  function  can  be  written  in  terms  of  the  mixture  velocity 
u  as  follows: 


exp 


=  exp 


1  (£-iy)2~ 

2  RaT 

1(£-«)21 

2 


x  exp 


RaT 

1  ( u 

2 


exp 

Ua)^ 


{£-u)-(u-ua) 


RaT 


RaT 


(B.50) 
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With  the  procedure  described  in  (4;  5),  fa  ^  of  Eq.  (B.22b)  is  the  second  order  Taylor 
expansion  of  the  first  exponential  in  the  right-hand  side  of  the  above  equality.  Because  the 
velocity  difference  (u  —  ua)  is  considered  to  be  a  small  quantity  (49),  the  second  exponential 
in  the  right-hand  side  of  the  above  equality  can  be  approximated  by  its  first  order  Taylor 
expansion: 

exp  [/3(t(£  -u)-(u-  ua)}  «  1  +  -it),  (B.51) 

where  /3a  =  l/RaT.  Since  the  velocity  difference  (it  —  ua)  is  small,  the  third  exponen¬ 
tial  in  the  right-hand  side  which  depends  on  ( u  —  «„)2  can  be  neglected  in  a  first-order 
approximation.  It  should  be  noted  that  this  term  affects  the  temperature  difference  [cf. 
Eqs.  (B.8b)  and  (B.12b)],  and  must  be  included  if  thermal  diffusion  is  important.  For  the 
case  of  isothermal  mixing  considered  here  the  simplifications  yield  the  final  result  of  the 
equilibrium  distribution  function  for  given  by  Eq.  (B.22a)  after  the  velocity  space  £ 
is  properly  discretized  (4;  5). 

B.A.3  Chapman-Enskog  Analysis 

In  the  Chapman-Enskog  expansion,  we  introduce  a  small  parameter  e  (which  is  the  Knudsen 
number),  and 


/a=/a(0)+*/a(1)  +  ---, 

(B.52a) 

dt  —  tdto  +  e2^ti  +  •  •  •  » 

(B.52b) 

V  eV. 

(B.52c) 
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The  left  hand  side  of  Eq.  (B.18)  is  expanded  in  a  Taylor  series  in  St  up  to  2nd  order  first, 
and  then  is  substituted  with  the  expansions  of  e: 

fcdxi  +  ea5t,t  +  St)  —  faixi>t) 

=  e8tDaf°  +  ^eH2tDaDara+--- 
=  eStD^f*® 

+e2St  +  d£0)/£(1)  +  \stD^D^f<0^ 

+  ■■■ 

where  Eqs.  (B.52)  have  been  substituted,  and 

Da  =  dt  +  ea- V,  and  =  dto  +  eQ-V. 

The  first  few  equations  of  a  set  of  successive  equations  in  the  order  of  e  obtained  from  the 
lattice  Boltzmann  equation  (B.18)  axe: 

e°:  /q(0)  =  /£(eq)  [l  +  (6a . —2'(^  ~  tt)1  ,  (B.53a) 

e1:  StD^f^  =  --/£(1)  +  J?  -  StFa,  (B.53b) 

Ta 

e2:  5t  (dhf^+D^f^  +  ^StD^D^f^ 

=  -^"/a(2)-  (B.53c) 

Ttr 

In  the  Chapman-Enskog  analysis,  it  is  assumed  that  the  mutual  interaction  term  J^q  as 
well  as  the  forcing  term  Fa  are  of  the  first  order  in  e  (17;  18;  63;  64). 

The  first  few  moments  of  the  equilibrium  distribution  function  can  be  easily  com- 
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puted: 


Q 

Y^fa°'>ea=  PaUa, 
a 

^  j  fix'  ^  eaieaj  =  RffT po-Sij  d"  p^UiUj 

a 

UiUjU-(ua  —  u) 


(B.54a) 

(B.54b) 


+P<T  |  (Ua  -  u){Uj  +  (ua  -  u)jU 

RtjTpvSij  PffU{Ujy 

£/Q(  =  RaTpa&ijklUl 

a 

RfjT P(jkijki{vi(j  u'ji  -f-  • 


1- 


R*T 
u-{ua  -  u) 


(B.54c) 


R„T 


RffR P<J  (ftijUpk  d"  fijk'U'ari  d*  ) > 


where  eQt-  is  the  i-th  Cartesian  component  of  the  vector  eQ  and 


(B.54d) 


^ ijkl  fiijtifd  d*  SilcSjl  d~  fiilfijk’ 

In  the  second  and  third  order  moments  of  fa ^  in  Eq.  (B.54c)  and  (B.54d),  respectively, 
the  terms  involving  products  of  the  velocity  and  the  velocity  difference  are  neglected. 

With  the  substitution  of  Eq.  (B.53b),  Eq.  (B.53c)  becomes: 


St 


^i/a(0)+40)/a(1)  +  \d£XJ?  -  5tFa) 


r 

To- 


(B.55) 


Therefore  the  second  order  (in  e)  equation  for  the  mass  conservation  is 


dtlP*  = 


1  Pops 


L  td  p 


(«<7  -  «?) 


=  -^iV  -pda. 


(B.56) 


where  the  spatial  variation  of  the  external  force  is  neglected  because  it  is  canceled  out 
by  properly  setting  the  velocity  change  due  to  the  external  forcing  to  5u  =  ^a5t  in  the 
equilibrium  distribution  function  (63;  64). 
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By  construction  of  the  forcing  term  Fa  [given  by  Eq.  (B.19c)],  we  have 

^  ^  =  0.  (B.57) 

a 

Therefore  the  forcing  term  has  no  contribution  to  the  pressure  tensor.  However,  the  partial 
pressure  tensor  for  each  individual  species  is  affected  by  the  cross-interaction  term  J 
because 

r/?e  .e  .  - 

/  v  Ja  “ 

a  rD  P 

x  ( UiSuj  +  UjSui)  -  — ^ -uiuju  •  Su  ,  (B.58) 

.  -tLcr-t 

where  Su  =  ( ua  —  uq),  and  Su{  =  ( ua{  —  The  effect  due  to  J £*  is  of  higher  order  in 
u  and  ttcr ,  and  thus  can  be  neglected  in  the  pressure.  Consequently,  the  partial  pressure 
tensor  can  be  approximated  as  the  following: 

a  a 

=  -rJt  dt0  fa{0)eaieaj  +  ^V- eaf^eaieaj 

,  a  a 

TpSt  [dt0  (RaT  Pc  Si  j  paUiUj )  R&T  {djPaUgi  +  dipa^aj  "4“  Sij1^  •  pcr'Uc r)]  +  •  •  • 

=  -Ta6t  [RaT(dt0Pa  +  'V'PcrUa)5ij  +  dto{paUiUj )  +  RaT  {djpauai  +  dipauaj)}  -I - 

~  -TffSfRaT  (djPvUvi  +  diPaUaj )  ,  (B.59) 

where  the  terms  smaller  than  0(M2)  (M  is  the  Mach  number)  are  dropped  as  usual  (65). 
Therefore, 

djPij1^  «  -TaStRaTpa^Ucri-  (B.60) 
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Appendix  C 


Computational  Assessment  of 
Diffusivity  in  LBM  for  Binary 
Scalar  Mixing 


Cs 

Speed  of  sound 

Nomenclature 

Cq! 

Discretized  particle  velocity 

c 

Lattice  constant 

D 

Diffusion  coefficient 

F 

Interaction  force 

Gf 

Interaction  strength  coefficient 

na 

Number  density  distribution 

Equilibrium  distribution  function 

m 

function 

Mass  of  molecule 

n 

Number  density 

n* 

Number  density  scalar 

p 

Pressure 

u 

Fluid  velocity 

p 

Fluid  density 

St 

Time  step 

sx 

Grid  size 

wa 

Weight  function 

T 

Normalized  relaxation  time 

C.l  Introduction 


In  Appendix  B,  a  LBE  model  for  binary  mixing  was  developed  directly  from  multi-component 
Boltmzann  Equation.  In  this  Appendix  we  will  develop  a  relatively  simple  model  by  adding 


83 


C.2.  Lattice  Boltzmann  model  for  binary  scalar  mixing 


84 


inter-particle  reaction.  This  model  has  a  simple  interaction  term  and  is  easy  to  be  imple¬ 
mented. 

The  remainder  of  this  Appendix  is  organized  as  follows.  In  section  C.2,  the  lattice 
Boltzmann  model  for  binary  scalar  mixing  is  presented.  Macroscopic  governing  equations 
for  each  species  are  derived.  Numerical  evaluations  are  presented  in  Section  C.3.  Finally 
we  conclude  in  Section  C.4  with  a  brief  discussion. 


C.2  Lattice  Boltzmann  model  for  binary  scalar  mixing 

In  a  two-species  field,  The  evolution  of  number  density  function  n£(£,  t)  (e.g.  a  =  1, 2)  for 
each  species  can  be  described  by  the  lattice  Boltzmann  equation 

K(x  +  eaSt,t  +  St)  -  =  — \[naa(x,t)  -  <(e?) (£,*)]  (C.l) 

T 

where  ra  is  the  relaxation  time  of  species  a;  {ea;a  =  0, 1,  •  ■  •  ,6}  is  the  discrete  velocity 
set.  For  the  sake  of  concreteness  and  simplicity  without  losing  generality,  we  shall  restrict 
our  derivation  to  a  nine- velocity  model  on  a  two-dimensional  square  lattice  (D2Q9  model). 
We  choose  the  equilibrium  distribution  functions  as 


rfW  =  NafaUrM)  =  wan[  1  + 


ua M)  9(ea 

O  "f” 


tp 


(e<?))2  3  (ua^f 


2c4 


2c2 


(C.2) 


where  ua(ei>')  is  a  equilibrium  velocity  which  is  determined  by  the  distribution  functions 
and  interaction  between  two  species.  The  sound  speed  of  this  model  is  cs  =  c/V 3,  where 
c  =  Sx/St  =  1  is  the  lattice  constant.  The  discrete  particle  velocities  ea  and  the  weighting 
factors  toQ  of  this  model  are 

\  0,  a  =  0 

(cos  ( )  c>  s{n  ( a  —  i_4 
(coS(i°L^)l  +  |)C)  +  |)c,  a  —  5_8 
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and 

4/9,  a  =  0 

wa  =  {  1/9,  a  =  1-4 

1/36,  a  =  5-8 

The  number  density,  mass  density,  velocity  of  the  cr-species  and  the  number  density,  the 

mass  density  and  velocity  of  the  mixture  are  calculated  by 


n<T =xx 

Q 

(C.3) 

pa  =  mcrna  —  rrf  ^  naa 

a 

(C.4) 

a 

(C.5) 

n  =  Y,n<7 

a 

(C.6) 

■Q 

II 

(C.7) 

u  =  -Y  pav? 

c 

(C.8) 

respectively. 

In  the  absence  of  interaction  forces  between  species,  we 

assume  that  the  velocities  u<7(e?) 

of  two  species  are  all  equal  to  a  common  velocity  u'.  Since  in  this  case  the  total  momentum 


of  particles  of  all  species  should  be  conserved  by  collision  operator  at  each  lattice  site,  it 
follows  directly  from  Eq  (C.l)  that 


C.2.  Lattice  Boltzmann  model  for  binary  scalar  mixing 


86 


Next,  we  introduce  interaction  force  between  two  species 

j*  =  <$nff/[Vx'  +  (x'-w')y]  (C.9) 

Where  Gf  is  a  factor  which  determines  the  strength  of  interaction;  x°  —  n<T/ n  is  number 
density  ratio;  and  oja  =  pa  /  p  is  mass  density  ratio.  To  incorporate  this  momentum  change 
in  the  dynamics  of  the  distribution  functions,  one  can  simply  define 

=  p°v!  +  Fa  (C.10) 

By  summing  Eq.(C.l)  over  all  directions  and  multiplying  Eq.(C.l)  by  ea  and  ma  and 
summing  over  all  the  nine  directions  and  two  species,  we  can  obtain  a  mass-conservation 
relation  for  each  species  and  momentum-conservation  relation  for  the  mixture  respectively. 
These  relations  are 


E™*^  +  eQ,t  +  1)  -  =  0  (C.ll) 

a  a 


E m<r  E * + 1)  -  e  E  <(*,  *)  =  E  = 0  (°-12) 

a  a  a  a  a 

C.2.1  Macroscopic  equations  of  each  species 

Following  the  algorithm  by  Xianwen  Shan  and  Gary  Doolen  (5;  6),  we  define  the  macroscopic 
velocity  of  the  whole  fluid  u  as  the  average  of  the  momentum  values  of  species  a  before  and 
after  collision  summed  over  two  species.  The  momentum  of  species  a  before  collision 


m 


(C.13) 


while  the  momentum  after  collision 

(1  -  X  ^  X  <We"«  (C.14) 
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so  that 

pu=l  Yip'tf  +  p'v?  -  =  pavF  +  \Fa  (C.15) 

2  1 f  T  T  2 

<J 

by  noticing  Eq.(C.lO). 

In  the  Chapman-Enskog  expansion,  we  introduce  a  small  parameter  e(which  is  the  Knud- 
sen  number),  and 

na  —  na(0)  +  ena(1)  H - >  dt  =  edto  +  e2dti  +  •  •  •  V  — >  eV  (C.16) 


Expanding  the  left-hand  side  of  Eq.(C.l)  in  a  Taylor  series  up  to  the  2nd  order  in  time,  we 
get 


n  t-7\  cr  .  1/q  ,  -•  y7\2  cr  _  [<(£,  t)  -  (n^  (S,t)] 

[dt  +  ea  •  V)na  +  -(dt  +  ea  •  V)  na  = - j -  (C.17) 

2  T 

Since  now  the  collision  operator  depends  upon  the  spatial  derivatives  of  particle  number, 
we  chose  the  leading  order  distribution  functions  to  be  the  equilibrium  distribution  about 
the  fluid  velocity  u,  namely,  n"Y  =  Na  (n17  ,u).  to  obtain  a  set  of  macroscopic  equations  in 
terms  of  the  correct  fluid  variables  without  changing  the  conservation  relations.  The  lattice 
Boltzmann  equation  is  then  satisfied  at  the  next  order  when  terms  that  depend  on  spatial 
derivatives  are  included.  Since 


E<{0)  = n<r’  Em<rEnS(0)r«  =  p*  (c-18) 

a  <7  ol 

we  get 


Cl  <t  ol  a 

Expanding  the  left-hand  side  of  Eq.(C.17)  by  (C.16),  we  have 
n°(x  +  ea,t  +  1)  -n°(x,t) 

~  e[(dtQ  +  ea  •  V)<(°)]  +  e2[(dt0  +  ea  ■  V)n^  +  8n <(0)  +  +  ea  •  V)2<(°)]  (C.20) 
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Substituting  above  equation  into  the  conversation  relations  (C.ll)  and  (C.12)  and  col¬ 
lecting  the  first  order  of  e,  the  following  equations  are  obtained 

dtopa  +  V  •  (pau)  =  0,  cr  =  1, 2  (C.21) 

dtQ(pu)  +  e^Vp+V  ■puu  =  YJFCT  (C.22) 

a 

Using  Eqs.(C.21),  we  rewrite  Eq.(C.22)  in  the  form 

{dto  +u-V)u  =  -ivP  (C.23) 

where  P  is  the  pressure  given  by 

P  =  C2p  +  £p-  (C.24) 

cr 

The  terms  of  the  2nd  order  of  e  in  the  expansion  of  the  mass-conversation  relation  yield  the 
following  equation  for  each  species 

dtiP*  +  mffV  •  Y,  <(1)?a  +  V  •  (dto  Y  <(0)e«  +  V-£  <(0)eaeQ)  =  0  (C.25) 

a  a  a 

Noticing  Eq.(C.21),  we  can  evaluate  the  third  term  of  Eq.(C.25) 

rrf(dto  Y  na  ' D)eQ  +  V  ■  Y  =  -w^Vp  +  c]Vpa  (C.26) 

a  a 

To  calculate  the  2nd  term  of  Eq.(C.25),  we  substitute  the  expansions  of  Eqs.(C.16)  into  the 
lattice  Boltzmann  equation  (C.l)  and  obtain 

dto  Y  na  +  V  •  Y nS( 0)eQea  =  --p[pa(u  -  v!)  +  Y  m^n^ea]  +  Fc  (C.27) 

cl  a  a 


which  when  multiplied  by  r0"  and  summed  over  all  the  components,  becomes 
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P(u  -  tf)  =  £  r  °f° + lYl  + E r<vr  -  °2s  E r<r  (c-28) 

o'  a  a  cr 

'^2m<Tn<a1^ea  can  then  be  solved  by  substituting  Eqs.  (C.28)  and  (C.26)  into  Eq.  (C.27). 

a 

Combining  Eq.  (C.25)  with  the  first  order  equation  (C.21),  we  obtain  the  following  equation 
at  the  second  order 

^  +  V  •  {p°u)  =  -r'V  •  Fa  +  (t°  -  i) V  •  [c2sVp°  -  w'Vp]  +  V  •  +  \)Fk 

k 

+' Vp  Tkuk  -  cl  rkVpk}  (C.29) 

k  k 

When  summed  over  two  species,  the  equation  above  becomes  the  continuity  equation  of  the 
whole  fluid  at  the  2nd  order 

^  +  v  •  (pu)  =  0  (C.30) 

C.2.2  Mutual  diffusivity  calculation 

The  species  of  a  fluid  mixture  will  be  diffused  into  each  other  if  the  mean  velocities  of  the 
species  differ.  The  local  velocity  of  the  fluid  mixture  can  be  defined  in  several  different  ways: 
by  averaging  the  velocities  of  the  constituents  by  mass,  mole,  or  volume(ll).  The  diffusion 
velocities  are  then  defined  relative  to  this  local  velocity.  Mathematically,  all  the  averaging 
methods  are  equally  useful  in  describing  the  diffusion  of  the  constituents.  Following  the 
treatment  presented  by  Chapman  and  Cowling(12),  we  use  the  mass  fluxes  of  the  species  in 
our  calculation.  Here  again,  since  the  momentum  of  each  species  changes  at  each  collision, 
we  must  average  mass  fluxes  before  and  after  collisions  to  get  the  overall  mass  flux 

P**’  =  ^l£<e„  +  (1  -  i)  (C.31) 

a  a  a 

By  applying  the  same  Chapman-Enskog  technique  above,  after  straightforward  manipulations^; 
6),  the  relative  mass  flux  of  species  a  is  obtained 

-  «")  =  -r’F’  +  (t"  -  |)[c *Vp°  -  w'Vy]  +  <^(£(7-*  + 

k 
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+Vp  £  Tkuk  -clY,  TfcVp*]  (C.32) 

k  k 

The  relative  velocity  n0"  —  is  the  mass-averaged  diffusion  velocity  of  species  a  with  respect 
to  u,  which  indicates  the  motion  of  species  a.  Noticing  that  the  right-hand  side  terms  of 
Eq.(C.29)  are  exactly  what  the  divergence  operator  acts  on  in  the  right-hand  side  of  above 
equation,  we  can  simply  rewrite  Eq.(C.29)  as  the  continuity  equation  of  species  a 

2£  +  V.(p'0  =  0  (C.33) 

which  demonstrates  that  each  species  satisfies  its  own  continuity  equation  at  2nd  order. 
Following  the  convention  in  the  diffusion  literature,  we  define  the  mass  flux  of  species  a  as 
ja  =  pa(u  —  ua).  As  a  result,  Eq.(C.29)  can  be  simply  expressed  as 

(dt  +  u-  V)p£r  +  V-u  =  -V-  f  (C.34) 

By  assuming  the  incompressibility  of  the  fluid(i.e.  V-it  =  0)  and  setting  same  molecular 
mass  ( ma  =  m)and  relaxation  time  (ra  =  r),we  obtain  the  following  advection-diffusion 
equation 

(dt  +  u-  V)p<r  =  -V  •  c2sp(rnGf  -  r  +  (C.35) 

The  diffusion  coefficient  in  the  mixture  is 

D  =  -c]p(TnGf  -  t  +  (C.36) 

C.3  Non-premixed  binary  scalar  mixing 

In  all  simulations  in  this  Appendix,  we  suppose  that  two  species  have  the  same  relaxation 
time  for  simplicity  and  periodic  boundary  conditions  are  used  in  both  directions. 

We  first  compute  the  mutual  diffusivity  in  a  binary  mixture  by  measuring  the  decay  of 
a  sinusoidal  concentration  wave  with  small  amplitude.  The  simulation  was  carried  out  on 
a  128  x  3  rectangular  lattice.  Initially,  we  assume  that  both  species  have  a  same  particle 
density  nb  =  nw  =  0.5  at  each  grid  in  the  whole  field.  A  perturbation  particle  density  wave 
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n(x)  —  0.01sin(2wx/128)  is  introduced  into  black  species,  i.e.  nb(x,  0)  =  0.5  +  n(x).  Two 
runs  with  r  =  1.0  and  r  =  2.0  are  performed. 

Analytical  solutions  are  calculated  from  Eq.  (C.36).  Fig.  C.l  exhibits  that  numerical 
computations  (dots)  are  in  excellent  agreement  with  the  analytical  predictions  (lines). 

Secondly,  we  study  the  binary  scalar  mixing  in  a  square  lattice.  Two  species  labeled 
black(b)  and  white(w)  are  initially  segregated  and  randomly  distributed  in  the  computa¬ 
tional  domain  as  shown  in  Figure  C.2. 

Initially,  the  macroscopic  velocity  is  set  to  zero  everywhere  in  the  domain  corresponding 
to  a  pure  diffusion  problem.  The  particle  densities  of  both  species  are  nb  =  1.0,  nw  =  0.0 
in  the  region  of  the  black  species  and  nw  =  1.0,  nb  =  0.0  in  the  region  of  white  species. 

It  should  be  reiterated  here  that  this  model  permits  the  simulation  of  scaler  miying 
between  species  of  equal  or  unequal  mass  molecular  weights.  The  mixing  characteristics  of 
both  cases  have  been  reported  in  the  paper  (10).  From  the  time  evolution  of  the  probability 
density  function(PDF)  of  number  density  scalar  n*  =  we  have  captured  the  changes 

of  the  PDF  shape  from  the  initial  double-delta  shape  to,  finally,  a  Gaussian-like  distribu¬ 
tion,  which  is  in  excellent  qualitative  agreement  with  DNS  data.  Here  we  further  test  the 
numerical  accuracy  of  the  model  for  equal-mass  case.  We  focus  on  comparisons  of  the  time 
evolutions  of  PDF  of  scalar  n*  with  different  interaction  strength  (diffusivity)  and  different 
resolutions  to  see  how  the  mixing  is  affected  by  these  parameters.  Three  snapshots  of  PDF 
shape  are  taken  in  each  run. 

Fig.  C.3  shows  the  time  evolutions  of  PDF  of  n*  with  different  interaction  strengths 
in  1282,  2562  and  5122  respectively.  Lines  correspond  to  Gf  =  0.002(D  =  0.166)  and 
symbols  to  Gf  =  0.125(D  =  0.125).  Three  time  instants  are  picked  up:  initial  double¬ 
delta  shape(dashed-dot  line  and  {>),  intermediate  uniform  distribution  (solid  line  and  A), 
and  final  Gaussian-like  shape  (dashed  line  and  *).  Time  instants  corresponding  to  initial, 
intermediate  and  final  have  been  normalized  by  a  dimensionless  constant  s  =  Dt/N  where 
N  denotes  the  length  of  the  computation  domain.  In  each  resolution,  PDF  shape  match 
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very  well  at  each  corresponding  mixing  stage,  which  demonstrates  that  the  mixing  process 
is  well  captured.  From  the  resolution  perspective,  5122  grid  performs  the  best  as  expected. 

C.4  Conclusions 

We  have  presented  a  LBM  model  for  binary  scalar  mixing  flows.  Macroscopic  equations 
governing  the  flows  are  derived  in  details  by  the  Chapman-Enskog  technique. 

Simulations  of  the  mixing  model  yield  encouraging  results:  (i).The  well  known  results  of 
PDF  evolution  from  continuum  Navier-Stokes  simulation  in  equal  mass  case  are  reproduced 
and  unequal  mass  case  generated  right  mixing  physics  which  is  difficult  to  handle  with 
continuum  based  numerical  method,  (presented  in  our  previous  paper(lO));  (ii).Mutual  dif- 
fusivity  are  calculated  analytically  and  confirmed  by  simulations;  (iii). Mixing  characteristics 
are  well  captured  in  different  resolutions. 

These  simulations  demonstrate  that  the  presented  LBM  models  are  practical  and  reli¬ 
able.  From  this  work,  more  complicated  flow  simulation  for  combustion  can  be  expected. 
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Figure  C.l:  Mutual  diffusivity  D  vs.  interaction  strength  Gf  for  given  rs  (r  =  1.0  upper  and 
t  =  2.0  down).  Lines:  analytical  prediction  from  Eq.(C.36);  Symbols:  LBM  computation. 


Figure  C.2:  A  typical  initial  number  density  distribution  of  two  species  colored 
with  black  and  white.  nb  =  1.0,  nw  =  0.0  in  region  of  the  black  species  and 
nw  =  1.0,  nb  =  0.0  in  region  of  white  species.  The  average  number  density  of 
both  is  nearly  the  same 


PDF 


Figure  C.3:  Evolutions  of  the  scalar  PDF  of  n*.  (a).  1282,  (b).  2562,  (c).  5122.  Lines: 
Gf  =  0.002(.D  =  0.166)  and  symbols:  Gj  =  0.125 (D  =  0.125).  Three  normalized  time 
stages:  initial(dashed-dot  lines),  intermediate(solid  lines)  and  final(dashed  lines) 
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Appendix  D 


Scalar  Mixing  and  Chemical 
Reaction  Simulations 


J 

Nomenclature 

Collision  operator  between  two  ea 

Discretized  particle  velocity 

na 

spices 

Number  density  distribution 

n(e?) 

Equilibrium  distribution  function 

m 

function 

Mass  of  molecule 

n 

Number  density 

u 

Fluid  velocity 

p 

Fluid  density 

St 

Time  step 

sx 

Grid  size 

wa 

Weight  function 

n 

Collision  operator 

xa 

Molar  fraction  for  spice  a 

u" 

Mass  fraction  for  spice  a 

C 

Concentration 

&OV 

Reaction  coefficient 

E 

Effective  activation  energy 

M 

Mass  weight 

R 

Universal  gas  constant 

Sl 

Burning  velocity 

T 

Temperature 

Y 

Mass  ratio 

D.l  Introduction 


With  few  notable  exceptions,  the  LBM  has  been  so  far  used  mostly  for  single-component,  in¬ 
ert  and  isothermal  flows.  In  this  Appendix,  we  simulate  scalar  mixing  in  a  multi-component 
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flow  and  a  chemical  reacting  flow  using  the  lattice  Boltzman  computational  approach.  At 
the  continuum  level,  the  mixing  example  considered  appears  as  a  pure  diffusion  problem 
without  any  advection  velocities.  However,  at  the  mesoscopic  level,  each  of  the  components 
is  associated  with  non-zero  velocities.  Hence,  the  problem  considered  is  truly  more  signifi¬ 
cant  than  simple  passive  scalar  mixing.  Due  to  the  kinetic  nature  of  the  LBE  scheme,  the 
extension  of  this  method  to  cases  with  non-trivial  macroscopic  velocity  field  is  straight  for¬ 
ward.  The  second  problem  considered  is  one-dimensional  flame  propagation  in  a  well-stirred 
homogeneous  mixture  of  propane  and  air.  The  flame  speed  calculated  using  the  LBM  is  in 
good  agreement  with  experimental  data.  This  problem  is  very  similar  to  the  one  solved  by 
Yamamoto(4)  but  with  a  sightly  different  physical  field.  The  overall  results  are  encouraging 
and  establish  the  feasibility  of  the  LBM  as  a  viable  computational  approach  for  simulating 
mixing  and  reaction. 

The  remainder  of  the  Appendix  is  organized  as  follows.  In  Section  D.2,  the  lattice 
Boltzman  equations  used  in  the  simulations  are  presented.  The  results  from  the  mixing 
and  reacting  simulations  are  presented  in  Section  D.4.  We  conclude  in  Section  D.5  with  a 
brief  discussion.  Details  of  the  physical  parameters  used  in  the  reacting  flow  simulation  are 
presented  in  Appendix  D.A.l 


D.2  Lattice  Boltzmann  equation  for  scalar  mixing 

For  the  sake  of  simplicity  without  losing  generality,  we  adopt  D2Q9  model,  as  shown  in 
Figure  A.l.  Consider  a  multi-species  field:  denotes  the  number  density  distribution 

function  of  a  particular  species  a  with  discrete  velocity  ea,  a  =  1,2,3,  .  The  number 

density  and  molecular  weight  of  species  a  are  given  respectively  by  na  and  ma .  Then  mass 
density  of  the  cr-species  is  given  as 

f7  ___  <7  \  A  __  (T 

p  —  771/  7%  —  771/  y  ^ 


(D.l) 
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Number  density,  mass  density  and  velocity  of  the  mixture  are 


n  =  Yjna 

a 

(D.2) 

p=Y,p° 

(D.3) 

and 

ff= 

P  a  a 

(D.4) 

respectively. 

For  the  sake  of  concreteness  and  simplicity  without  losing  generality,  we  construct  a 
LBE  for  binary  mixture  with  6X  —  St  —  1.  Let  a  and  a'  denote  the  two  species  of  interest. 


The  LBE  for  each  species  is 


n°{x  +  ea,t  +  1)  =n°{x,t) +  £l°{x,t) 
where  the  collision  operator 

1  jaaf 

ns  =  ~K  -»?•*>]  +  J-*~ 

ra  ma 

includes  an  additional  term  J which  reflects  interactions  between  two  species. 
We  use  the  following  number  density  equilibrium  distribution  function 


(D.5) 


(D.6) 


=  warf[\  +  3(ea  ■  if)  +  ~(ea  •  u')2  -  -u'2} 


(D.7) 


with 


u 


p 


The  binary  interaction  term  is  modeled  as  (3) 

1 


J° 


T^a1  •  [Va;*7  +  (xa  —  w*7)— ]  =  —  Jq,<t 


TUU  -  -  p 

where  xa  and  u>a  are  molar  and  mass  fractions  of  the  species  a 


xa  = 


rf 

n 


ua  =  ?- 


and 


n(ag)  =  wan[l  +  3(ea  •  u)  +  |(ea  •  u)2  -  ^u2] 


(D.8) 


(D.9) 


(D.10) 
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D.3  Lattice  Boltzmann  equations  for  reacting  flow 

Here  we  simulate  simple  one-dimensional  flame  propagation  through  a  premixed  mixture  of 
propane  and  air.  The  problem  studied  is  identical  to  that  of  Yamamoto(4)  but  the  physical 
combustion  field  is  sightly  different.  The  simplifying  assumptions  invoked  in  this  study  are 
now  listed: 

•  No  external  forces  in  the  field  and  the  flow  is  incompressible. 

•  The  chemical  reaction  (heat  release)  does  not  affect  the  flow  field;  temperature  and 
concentration  fields  are  decoupled  and  solved  separately. 

•  Nitrogen  is  inert. 

•  The  transport  properties  are  constants  (not  functions  of  temperature). 

•  Viscous  energy  dissipation  and  radiative  heat  losses  are  negligible. 

•  Simple  one  step  reaction  is  considered 

-f-  5O2  3(702  d"  4//2 O 

and  the  over-all  reaction  rate  is  given  by 

Uov  =  KovCc3HsCo2e~E/RT 

where  Cc3h8i  Co2>  Kov  and  E  are  concentrations  of  fuel  propane  and  oxygen,  reaction 
coefficient  and  effective  activation  energy  respectively. 

In  a  reacting  flow,  the  state  of  the  fluid  at  any  given  point  in  space  and  time  can 
be  completely  specified  in  terms  of  fluid  velocity,  composition  vector  (either  in  terms  of 
mass  fraction  or  concentration)  and  temperature.  We  will  need  to  develop  the  LBE  for 
all  these  variables.  For  generating  a  background  flow,  the  conventional  LBM  substeps  of 
collision(relaxation)  and  streaming  (convection)  are  used.  However  for  the  temperature  and 
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concentration  fields,  there  is  an  extra  substep  between  collision  and  streaming  substeps  to 
account  for  reaction.  This  is  identical  to  the  time-splitting  approach  used  in  continuum 
methods  for  chemically  reacting  flows. 

Flow  field.  The  background  flow-field  is  obtained  using  the  following  stencil  for  partial 
pressure 

pa{x  +  ea,t  +  l)  =pa(x,t)  -  —\pa(x,t)  -pe^(x,t)]  (D.ll) 

Tp 

where 

Pa  =  WM 1  +  3(ea  •  u)  +  ^(ea  ■  uf  -  ^ u 2]  (D.12) 

The  total  pressure  p(=  pc2)  and  the  fluid  velocity  are  calculated  using 

P  =  J^Pa 

a 

U  =  “  Y]  ZaPa 

T)  Z ' 
r  Q 

This  is  the  velocity  used  for  determining  the  equilibrium  distribution  functions  in  temper¬ 
ature  and  concentration  fields. 

Temperature  and  concentration  fields.  For  temperature  (T  is  normalized  by  Tc  = 
E/R)  and  concentration  (mass  ratio  Y\i  =  CsHs,02,C02  and  H20)  fields,  there  is  an 
extra  computational  substep,  reaction,  besides  conventional  computational  substeps  of  col¬ 
lision  and  advection. 


•  COLLISION 


Ta(x,t)  =  Ta(x,t)  -  ~[Ta(x,t)-T^(x,t)} 

TT 

(D.13) 

Km = -  —  K(*,<)  -^<e3)(s,<)] 

Tyi 

(D.14) 

?ieq)  =  waT[  1  +  3(ea  ■  u)  +  \{ea  ■  uf  -  ^u2} 

(D.15) 

Yla{eq)  =  Wa^i  1  +  3(eQ  •  u)  +  \(ea  ■  uf  -  ^u2] 

2  2 

(D.16) 

where 
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and 

a  a 

Relaxation  time-constant  tt  is  determined  by  thermal  diffusivity  and  rY% ’s  are  deter¬ 
mined  by  the  diffusivity  of  corresponding  species. 


REACTION 
Reaction  equation 

CSHS  +  5 02  ->  3C02  +  4 H20 
a^ov  K>ovCc$H%C'02i^  ' 

Concentrations 

Ci  =  pYi/Mi 


Reaction  terms 


where 


r\  _  \  o 

P 

Vt  =  — 7pruov 
pCpTc 


**  =  E%  f=Ef« 

a  a 

In  the  above  equations,  the  stoichiometric  coefficients  (A’s)  for  the  various  species  are: 
Ac3/f8  =  —  I,  A o2  =  —5,  Aco2  =  3,  A#2o  =  4  and  all  physical  parameters  are  listed  in 
Appendix  A. 


•  STREAMING 


Ta(x  +  ea,t  +  1)  =  Ta(x,t)  +  waQT 
Y^{x  +  ea,  t  +  1)  =  y^(f,  t)  +  waQYi 
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D.4  Simulations 

As  mentioned  in  the  introduction,  the  primary  objective  of  this  work  is  to  investigate  the 
ability  of  the  LBM  to  simulate  scalar  mixing,  chemical  reaction  and  ultimately  turbulent 
combustion.  Working  towards  this  end,  we  perform  simulations  of  two  unit  problems:  one 
to  establish  the  mixing  characteristics  and  another  to  demonstrate  the  chemical  reaction 
scheme. 


D.4.1  Non-premixed  binary  scalar  mixing 

This  problem  epitomizes  the  scalar  mixing  issues  encountered  in  a  typical  non-premixed 
combustion  application.  Two  species  (presumably  fuel  and  oxidizer)  are  initially  segregated 
and  randomly  distributed  in  the  computational  domain  which  in  the  present  case  is  a  square 
box.  Mesh  size  is  set  500  x  500.  The  two  species  are  generically  labeled  as  black  and  white. 
A  typical  initial  distribution  is  shown  in  Figure  D.l. 

The  macroscopic  velocity  is  set  everywhere  to  zero  corresponding  to  a  pure  diffusion 
problem.  It  should  be  reiterated  here  that  the  mesoscopic  velocities  are  non-zero.  The 
initial  values  for  the  number  densities  are  nb  =  1.0,  nw  =  0.0  in  region  of  the  black  species 
and  nw  —  1.0,  nb  =  0.0  in  region  of  white  species.  Simulations  are  performed  for  two 
different  sets  of  relaxation  time  scales:  r6  =  tw  =  1.0  and  rbw  =  rwb  =  30.0.  Citing 
homogeneity  of  the  scalar  field,  periodic  boundary  conditions  are  used  in  all  directions. 

To  discretize  Vxa  and  Vp  in  the  binary  interaction  term  in  Equation  D.9,  we  use  central 
difference  operator.  Taylor  series  expansion  of  each  of  /(£+$)  terms  to  the  second  order 
leads  the  following  stencil  of  dXl  and  dX2 : 


-1 

0 

1 

„  1 

1 

4 

1 

-4 

0 

4 

’  9x2  ~  12 

0 

0 

0 

\ - 

1 

l_i 

0 

1 

-1 

-4 

-1 

(D.17) 
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That  is 

8xi /(*)  ~  ^[4/(f+ei)  +  /(f+e5)  +  /(f+e8)-/(x  +  e6)-4/(s  +  e3)-/(f+e7)]  (D.18) 

dx2f(%)  ~  ^[4/ (a?+  e^)  +  f(x  +  £5)  +  f(x  +  e$)  —  f(x  +  eV)  —  4 /(x  +  64)  —  /(x  +  eg)]  (D.19) 
in  which  /  can  be  x^  and  pa. 

The  lattice  Boltzmann  methodology  permits  simulation  of  mixing  between  species  of 
equal  or  unequal  number  and  mass  densities  with  equal  facility.  However,  in  continuum 
based  methods,  mixing  between  species  of  unequal  mass  densities  is  not  straight-forward. 
This  represents  a  fundamental  advantage  of  the  LBM  over  continuum-based  methods. 
Equal  mass  case(m6  =  mw  =  1.0).  The  first  case  studied  is  mixing  of  two  fluids  of 
equal  molecular  weights  and  number  densities,  hence  of  equal  mass  density.  This  case  is 
particularly  interesting  as  the  results  can  be  directly  compared  with  directed  numerical 
simulation(DNS)  of  Navier-stokes  equation  data  of  Eswaran  and  Pope  (5).  In  this  case, 
the  number  density  and  mass  density  are  equivalent  since  the  molecular  weight  of  the  two 
species  are  identical. 

Figure  D.2(a)  shows  the  time  evolution  of  the  probability  density  function(pdf)  of  scalar 

p: 


The  corresponding  DNS(5)  data  is  shown  in  Figure  D.2(b).  The  LBE  and  DNS  data 
show  excellent  qualitative  agreement.  In  particular,  the  change  of  the  pdf  shape  from  the 
initial  double-delta  shape  through  a  nearly  uniform  distribution  to,  finally,  a  Gaussian-like 
distribution  is  well  captured  by  the  LBE  results.  It  deserves  mention  here  that  many  other 
mixing  models  do  not,  even  qualitatively,  capture  the  form  of  pdf  during  evolution. 

In  Figure  D.3,  the  time  evolution  of  the  root-mean-square(rms)  of  scalar  fluctuations 
(p1)  obtained  from  LBE  is  compared  with  that  from  DNS  (5): 


p'  =  \f  <  (pb~  <  pb  >)2  > 


(D.20) 
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where  <  •••  >  implies  volume-averaged  value.  The  relaxation  time-constant  has  been 

chosen  to  yield  the  best  agreement.  For  the  optimal  time-constant,  the  agreement  is  again 

excellent.  In  ongoing  work,  we  attempt  to  accurately  characterize  the  relationship  between 

the  time-constant  and  diffusion  coefficient  (3).  Here,  it  suffices  to  say  that  given  the  right 

time-constant,  the  LBE  captures  the  DNS  behavior  well,  qualitatively  and  quantitatively. 

Unequal  mass  case(m6  =  2.0,  mw  =  1.0).  The  unequal  mass  case  is  particularly 

interesting  since  it  represents  a  more  practical  problem,  mixing  between  species  of  unequal 

mass  densities.  In  this  case,  the  initial  distribution  of  black  and  white  species  is  similar 

to  the  equal  mass  case.  Thus  the  average  number  density  of  black  and  white  species  are 

identical.  However,  the  molecular  weight  of  the  black  species  is  twice  that  of  the  white 

species.  Hence,  the  macroscopic  mass  density  of  the  black  species  is  twice  that  of  the  white 

species.  The  precise  definition  of  the  mass  density  used  here  is 

—  _  pa  —  pa'  _  rrfrf  —  m^n*' 

^  pa  +  pa'  m?na  -I-  m?'  r\F' 

and  the  particle  number  density  is 

— -  na  -  na' 

na  = - T 

na  +  n° 

In  the  unequal  molecular  weight  case,  mass  density  and  number  density  are  two  independent 
entities.  The  evolution  of  both  these  quantities  are  investigated.  In  Figure  D.4,  the  mass 
density  evolution  is  given  for  both  species.  The  mass  density  goes  from  an  initial  double¬ 
delta  shape  to  a  Gaussian-like  shape  centered  around  the  global  average  of  the  respective 
density.  The  final  pdf  shape  for  each  species  is  clearly  not  symmetric.  This  is  easily 
understood  since  the  overall  average  density  of  the  black  species  is  twice  that  of  the  white 
species. 

In  Figure  D.5,  the  number  density  evolution  is  shown.  Since  the  average  number  density 
is  identical  for  both  species,  the  final  pdf  distribution  is  symmetric  about  the  mean  value 
for  each  species.  However,  the  intermediate  forms  of  the  pdf  are  quite  nonsymmetric, 
demonstrating  that  the  mixing  process  in  this  case  is  quite  different  from  the  equal-mass 
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case  even  if  the  final  distributions  are  similar.  A  detailed  study  of  the  physics  of  mi-ring 
between  species  of  unequal  mass-densities  will  be  undertaken  later. 

D.4.2  Reacting  flow  in  a  1-D  channel 

In  this  example,  we  study  the  ability  of  LBM  to  simulate  chemical  reaction.  The  simplest 
non-trivial  case  when  reaction  can  be  studied  without  the  complicating  effects  of  mixing  is 
the  case  of  1-D  flame  propagation  through  a  homogeneous  premixed  mixture. 

Schematic  of  the  flow  simulated  is  shown  in  Figure  D.6.  For  this  simple  case,  the 
background  flow  generated  from  Eq.(D.ll)  maintains  both  the  pressure  and  velocity  fields 
uniform  in  space  and  time.  A  heat  source  is  placed  at  a  location  close  to  the  inlet  to  ignite 
the  mixture.  Once  ignition  is  achieved,  the  heat  source  is  removed.  At  subsequent  times 
the  flame  propagates  to  the  right.  Initial  conditions  are  set  as  following: 

The  values  of  pressure  and  velocity  are  set  at  p  =  1,  ux  =  uin  =  0.1,  uy  =  0.0.  Both  fields 
are  maintained  uniform  at  all  times  in  this  simple  case.  The  temperature  is  set  at  T  =  300° K 
everywhere  except  at  x  =  50  where  a  heat  source  is  placed  with  Tsource  =  1500° K  to  ignite 
the  mixture.  The  hot  spot  is  removed  after  the  mixture  ignites.  The  mass  ratio  of  nitrogen 
is  Yjy2  =  0.7375.  The  well-premixed  mixture  consists  of  propane  and  oxygen  with  the  mass 
ratios  of  Yc3hs  —  0.2252,  Yo2  =  0.0373.  The  mass  fractions  of  the  products  are  initially  set 
to  zero:  YCo2  =  YHio  =  0.0. 

Periodic  boundary  conditions  are  used  at  the  top  and  bottom  boundaries  and  the  fully 
developed  boundary  condition  is  applied  at  the  outlet.  At  the  inlet,  the  initial  conditions 
are  maintained. 

In  Figure  D.7,  the  flame  position  is  shown  as  a  function  of  time.  The  flame  location 
is  identified  as  the  position  with  the  highest  reaction  rate  at  any  given  time.  The  linear 
variation  of  flame  location  with  time  (in  Figure  D.7)  indicates  that  the  flame  propagates  at 
a  nearly  constant  rate.  This  flame  speed  can  be  easily  estimated  from  knowing  the  flame 
position  at  initial(fj  =  0;  Xfi  =  50)  and  final(fy  =  4000;  Xfi  =  406)  times.  The  flame  speed 
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thus  calculated  is 


f~  tf-U 


406  -  50 
4000  -  0 


0.0089 


in  lattice  units.  Knowing  the  flame  speed,  the  burning  velocity  can  be  easily  determined: 


^ L  —  U{n  V  j 

In  the  above,  Uin  is  the  reactant  velocity  at  the  inlet(which  is  maintained  uniform  throughout 
the  flow-field).  The  burning  velocity  thus  obtained  will  be  in  lattice  units.  This  can  be 
converted  into  metric  units  as  follows: 


Sl  = 


Uin  Vf 
Uin 


Up 


(D.21) 


The  resulting  burning  velocity  is  SL  =  0.11  m/s  which  compares  extremely  well  with  the 
value  obtained  from  experiments  for  a  propane-air  flame(6). 

Figure  D.8  shows  that  the  reaction  rate  profile  in  the  reaction  zone  as  time  evolves. 
Simulations  indicate  that  flame  behavior  is  sensitive  to  the  magnitude  of  the  heat  source. 


D.5  Conclusions 

In  this  Appendix,  we  simulate  scalar  mixing  and  chemical  reacting  flows  using  LBM.  In 
the  case  of  equal-density  species  mixing,  well  known  results  from  continuum  Navier-stokes 
simulation  are  reproduced.  The  true  advantage  of  the  LBM  can  be  seen  from  the  mixing 
simulations  of  species  of  different  molecular  weights.  The  results  appear  quite  encouraging. 
Such  simulations  are  very  difficult  with  continuum  based  methods.  The  premixed  reacting 
flow  simulations  also  produce  results  that  are  in  good  agreement  with  known  data.  Based 
on  these  simulations,  we  conclude  that  LBM  can  perform  adequately  for  more  complicated 
turbulent  combustion  simulations. 


Figure  D.4:  Pdf  evolution  of  mass  density  in  unequal  mass  case(m 


D.5.  Conclusions 


Figure  D.5:  Pdf  evolution  of 


=Uo 


D.6.  Figures 


112 


D.6  Figures 

Appendix  D.A.l:  Physical  parameters  used  in  the  reacting 
flow  simulation 


•  Reaction  coefficient  kov  =  9.9  x  107[m3  ■  mol -1  ■  s-1] 

•  Universal  gas  constant  R  =  8.315[J  •  mol -1  •  K-1] 

•  Effective  activation  energy  E  =  30 [kcal  ■  mol-1]  =  1.26  x  105[J  •  mol-1] 

•  Heat  of  overall  reaction  Q  =  2.05  x  106[J  •  mol-1] 

•  Density  of  pre-mixed  mixture  p  =  1.2\kg  •  m-3] 

•  Heat  capacity  Cp  =  29.1[J  •  mol-1  ■  K-1]  =  103[J  •  kg-1  ■  K-1] 

•  Kinetic  viscosity  v  =  1.6  x  10_5[m2  •  s_1] 

•  Thermal  diffusivity  k  =  2.2  x  10-5[m2  ■  s-1] 

•  Diffusivities 

Dc3Ha  =  1-1  x  10_5[m2  •  s-1],  Do2  =  2.1  x  10“5[m2  •  s-1] 
Dco2  =  1-6  x  10-5[m2  •  s-1],  Dh2o  =  2.2  x  10-5[m2  •  s-1] 

•  Mass  weight 


Mc3h8  =  4.4  x  10  2 [kg /mol],  Mo2  =  3.2  x  10 ~2[kg/mol] 
Mqo2  =  4.4  x  10 -2[kg/mol],  MHi0  =  1.8  x  10-2[kg/mol] 


•  Equivalent  ratio 


,  _  Yc3Ha/Yp2 
9  0.276 


=  0.6 


•  Length  of  the  1-D  channel  L  =  16.7[mm] 

•  Physical  velocity  up  =  1.0[m  •  s-1] 
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Appendix  E 


DNS  and  LES  of  Decaying 
Isotropic  Turbulence  With 
and  Without  Frame  Rotation 


Nomenclature 


Csm  Smagorinsky  constant 

cs  Speed  of  sound 

/  Distribution  function 

u  Fluid  velocity 

v  Viscosity 

St  Time  step 

r  Normalized  relaxation  time 

O,  Angular  velocity 

k  Wave  number 

£  Dissipation  rate 

Re\  Taylor  microscale  Reynolds  num¬ 
ber 


a  Coriolis  acceleration 

ea  Discretized  particle  velocity 

f(ei)  Equilibrium  distribution  function 

p  Fluid  density 

Sij  Strain  rate 

Sx  Grid  size 

wa  Weight  function 

n  Decay  exponential 

k  Turbulent  kinetic  energy 

A  Taylor  microscale  length 

t'  Normalized  time  ( teo/ko ) 
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E.l  Introduction 

In  this  Appendix,  we  explore  the  application  of  LBM  to  turbulence.  In  an  effort  to  evaluate 
the  capability  of  the  LBM  in  turbulence,  we  perform  direct  numerical  simulation  (DNS) 
and  large-eddy  simulation  (LES)  of  decaying  homogenous  isotopic  turbulence  (HIT)  in 
both  inertial  and  rotating  frames  of  reference.  Decaying  HIT  is  an  important  benchmark 
problem  in  the  field  of  DNS  and  LES  of  turbulence.  In  fact,  the  first  attempt  at  DNS 
with  incompressible  Navier-Stokes  equation  involved  this  problem  (1).  Since  then  several 
numerical  investigations  of  decaying  HIT  have  been  carried  out,  including  some  recent 
NS  studies  on  decay  exponents  and  low  wave-number  spectral  scaling  (2;  3;  4;  5).  Some 
preliminary  studies  of  three-dimensional  (3D)  decaying  HIT  using  LBE  have  also  been 
performed  (6;  7;  8),  but  these  investigations  stop  well  short  of  quantitative  comparisons 
with  the  well  established  classical  results. 

The  objective  of  the  present  Appendix  is  to  perform  a  comprehensive  investigation  of  de¬ 
caying  HIT  with  LBE-DNS  and  LBE-LES  to  establish  the  suitability  of  the  LBM  for  turbu¬ 
lence  applications.  For  this  purpose,  we  perform  three  type  of  simulations:  (a)LBE-DNS  of 
decaying  HIT  in  inertial  and  rotating  frame  of  reference.  The  decay  exponent  for 
the  kinetic  energy  k  and  the  dissipation  rate  e  are  computed  and  compared  with  correspond¬ 
ing  NS-DNS  results.  The  low  wave-number  scalings  of  the  energy  spectrum  are  studied.  The 
effect  of  rotation  on  the  kinetic  energy  decay  is  investigated,  (b)  LBE-LES  of  decaying 
HIT  inertial  frame  of  reference.  We  compute  kinetic  energy  decay,  energy  spectrum 
and  flow  structures  using  LBE-LES.  By  comparing  LBE-LES  results  with  the  corresponding 
LBE-DNS  results,  we  observe  that  LBE-LES  accurately  captures  large  scale  flow  behavior. 
We  find  that  the  optimal  Smagorinsky  constant  value  for  LBE-LES  is  smaller  than  the  tra¬ 
ditional  value  used  in  NS-LES  approaches,  (c)  LBE-LES  vs.  NS-LES.  We  carry  out  a 
comparative  study  of  the  LBE-LES  and  NS-LES  of  decaying  HIT.  We  show  that  the  LBE- 
LES  simulations  preserve  flow  structures  more  accurately  than  the  NS-LES  counterpart. 
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This  is  due  to  the  fact  that  some  history/non-local  effects  are  inherent  in  the  LBE  subgrid 
closure. 

The  remainder  of  this  Appendix  is  organized  as  follows.  Section  E.2  briefly  reviews 
relevant  background  on  decaying  HIT.  Section  E.3  gives  a  concise  introduction  to  the  LBM 
equations  for  DNS  and  LES.  We  present  our  results  in  Section  E.4  and  conclude  in  Sec¬ 
tion  E.5. 


E.2  Homogeneous  isotropic  turbulence 


The  energy  spectrum  E(k,  t)  in  decaying  HIT  evolves  as 


dtE(K,t )  =  — T(/c,t)  —  2  i>K2E(K,t), 


(E.l) 


where  k  is  the  wave-number  and  v  is  the  kinematic  viscosity,  and  T(k,  t)  represents  the 
nonlinear  energy  transfer  between  modes  (cf.  Eq.  (6.162)  in  (9)).  The  kinetic  energy  k  and 
dissipation  rate  e  of  turbulence  are  given,  respectively,  by 


k  =  J  E(k)(1k ,  and  e  =  2v  J  k2E(k)cIk. 

It  has  been  long  observed  that,  after  a  short  initial  transient  period  of  time,  the  kinetic 
energy  k  and  dissipation  rate  e  exhibit  power-law  decay  (9) 


where  ko  and  eo  are  the  values  of  k  and  e  at  the  reference  time  to  =  nko/eo.  Isotropic 
turbulence  is  typically  characterized  by  the  Taylor-microscale  Reynolds  number 

**»- *?* -!/!;*•  (E-3) 
where  A  =  ^15vu^ms/e  is  the  transverse  Taylor-microscale  length  and  urms  =  y/2k/3  is  the 
root  mean  square  (rms)  of  the  velocity  field  u. 

Equation  (E.l)  admits  a  continuous  class  of  invariant  solutions  in  the  limit  of  Re  — »•  oo 
(10).  At  the  large  Re,  E(K,t)  at  the  low  wave-number  behaves  as  limK_>o  E(k)  ~  k0',  where 
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cr  is  a  time-independent  constant  (e.g.  (11)).  For  inviscid  fluids,  if  Loitsyansky’s  integral 
(12)  is  an  invariant,  then  a  =  4  and  n  =  10/7  (13);  if  BirkhofF’s  integral  (14)  is  an  invariant, 
then  a  =  2  and  n  =  6/5  (15).  It  has  been  recently  shown  that  time-invariant  integral  length 
scale  l  corresponds  to  cr  =  oo  and  n  =  2  and  time-invariant  Reynolds  number  corresponds 
to  a  =  1  and  n  =  1(16).  Furthermore,  the  conservation  of  energy,  angular  momentum, 
and  helicity  lead  to  u  —  2,  7,  and  1,  in  the  limit  of  Re  — t  oo,  respectively.  The  energy 
conservation  of  inviscid  fluid  uniquely  determines  the  invariant  solution  of  Eq.  (E.l),  i.e., 
a  —  2,  in  accordance  with  Birkhoff’s  invariant  (14).  Despite  the  apparent  simplicity  of  the 
decaying  HIT  problem,  the  relevant  flow  invariant,  asymptotic  decay  exponent  and  the  low 
wave-number  scaling  are  strong  functions  of  the  initial  spectrum  and  Reynolds  number. 
There  is  still  no  clear  consensus  on  whether  the  angular  momentum  or  energy  is  the  correct 
invariant.  It  is  also  not  clear  what  the  conditions  are  under  which  the  invariance  of  either 
quantity  can  be  observed.  Consequently,  various  results  have  been  reported  (2;  3;  10;  17). 

We  perform  detailed  comparisons  with  established  data  qualitatively  and  quantitatively 
on  the  following  important  items:  (i)  energy  decay  exponent  n,  (ii)  low  wave-number  scaling 
of  the  spectra,  (iii)  flow  structure,  and  (iv)  effect  of  rotation  on  kinetic  energy  decay. 

E.3  LBE  formulation  for  DNS  and  LES  of  turbulence 

E.3.1  Lattice  Boltzmann  equation  for  DNS 

The  LBE  with  single-relaxation-time  approximation  due  to  Bhatnagar,  Gross,  and  Krook 
(BGK)  (18)  for  the  collision  operator  is  (19;  20) 

fct(x  +  ect^ii  t  +  St)  =  fa(x ,  t)  —  —  ^/Q  —  +  Fa,  (E.4) 

where  fa  is  the  density  distribution  function  with  discrete  velocity  ea  along  the  ath  direc- 
tion,  fa  is  the  equilibrium  distribution  function,  and  t  is  the  relaxation  time  due  to  the 
fluid  particle  collision  determining  the  viscosity  v  of  the  modeled  fluid.  In  what  follows,  we 
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use  the  LBE  model  with  19  velocities  in  three  dimensions,  i.e.,  the  D3Q19  model  shown  in 
Fig.  A.2. 


The  discrete  velocities  are: 

'  (0,  0), 

ea  =  ( 


(±1,  0,  0)c,  (0,  ±1,  0)c,  (0,  0,  ±l)c, 

(±1,  ±1,  0)c,  (±1,  0,  ±l)c,  (±1,  ±1,  0)c, 

The  equilibria  for  incompressible  flow  (21)  are 

[" 3ea  •  u  ^  9 (eQ  •  it)2  3 u2  "|  | 


a  =  0 
a  =  1-6 
a  =  7-18. 


(E.5) 


faq)  ~  W>a  | 


Sp  +  po 


(E.6) 


2c4  2c2  J  /  ’ 

where  dp  is  the  density  fluctuation,  and  po  is  the  constant  mean  density  in  the  system  which 
is  usually  set  to  1,  and  c  =  6x/5t  =  1  in  lattice  units  (i.e.  6t  =  Sx).  The  sound  speed  of 
the  model  is  cs  =  c/ \/3.  The  total  density  is  p  =  po  +  Sp.  The  weighting  factors  wa  for  the 
D3Q19  model  are  wo  =  1/3,  wi-q  =  1/18,  and  107-18  =  1/36.  The  mass  and  momentum 
conservations  are  strictly  enforced: 


OL  Ct 

(E.7a) 

P0U  =  ^2  eafa  =  e»faq)- 

ct  a 

(E.7b) 

For  athermal  fluids,  the  forcing  term  Fa  is  (22) 

Fa  =  -3wapoeQ'att, 

cz 

(E.8) 

where  a  is  the  acceleration  due  to  external  force.  In  our  simulations,  only  for  rotating  case, 

we  consider  the  Coriolis  force,  i.e.,  a  =  — 2f2  x  u,  where  Cl  is  the  angular  velocity  of  the 
frame  of  reference. 

The  hydrodynamic  equations  derived  from  Eq.  (E.4)  via  the  Chapman-Enskog  analysis 


are 


dtp  +  V  •  pu  =  0, 

dtu  +  u  •  V«  =  —  Vp  +  vV2u  +  a, 


(E.9a) 

(E.9b) 
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where  p  =  c2sp/p0  and  the  kinematic  viscosity  v  has  the  following  relation  with  the  relaxation 
time 


(E.10) 


It  is  important  to  note  that  in  LBE  the  strain  rate  tensor  S ij  can  be  obtained  directly 
from  the  second-order  moment  of  the  nonequilibrium  distribution  function  as 


S  ij  — 


2po 


)C2r  zL,  eaieaj  [/o  5 


so  that  the  dissipation  rate  is  computed  as  e  =  2u  .  SjjSjj. 


(E.ll) 


E.3.2  LES  extension  of  lattice  Boltzmann  equation 

For  the  LES,  high  wave  number  Fourier  components  of  the  density  distribution  function  are 
filtered  and  the  resolved-scale  distribution  function  is  separated  from  the  unresolved  part. 
The  filtered  form  of  the  LBE  for  LES  is  modeled  as  (23): 

7a(x  +  eaSt,  t  +  6t)  =  fa(x,t)  7ieq)]  +  Fa,  (E.12) 

where  / a  and  / Q  represent  the  distribution  function  and  the  equilibrium  function  of  the 
resolved  scales  respectively.  The  effect  of  the  unresolved  scale  motion  is  modeled  through 
an  effective  collision  which  has  been  included  in  the  LES  effective  relaxation  time  r*  in 
Eq.  (E.12).  The  LES  effective  viscosity  v *  is  then  obtained  from 

=  V  +  Vt  =  i  (t*  -  i)  cdx,  (E.13) 

with  vt  denoting  turbulent  viscosity  usually  called  eddy  viscosity. 

To  evaluate  the  fidelity  of  the  LBE-LES  simulations,  we  use  the  Smagorinsky  model 
(9;  24)  for  the  small  unresolved  scale  motion.  In  the  Smagorinsky  model,  the  eddy  viscosity 
vt  is  calculated  from  the  filtered  strain  rate  tensor  S tj  =  ( djUi  +  diUj)/ 2  and  a  filter  length 
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scale  Sx  as  follows: 

vt  =  (C'sm^x)2S,  (E.14a) 

(E.14b) 

where  S  is  the  characteristic  filtered  rate  of  strain  and  Csm  is  the  Smagorinsky  constant. 
With  Csm  and  6X  given,  Tt  can  be  obtained  from  Eq  (E.ll)  (23): 

Tt  =  \  (V t2  +  18V2(p0c2)-^C^mSxS  -  .  (E.15) 

As  shown  in  Eq.  (E.ll),  the  filtered  strain  rate  tensor  S ij  can  be  computed  from  the  second- 
order  moment  of  the  filtered  nonequilibrium  distribution  function  directly.  In  Eq.  (E.15), 
the  filtered  strain  rate  S  which  is  used  to  determine  vt  is  current  in  time.  Because  the  time 
step  in  the  LBE  simulations  is  relatively  small  in  physical  units,  we  can  also  use  the  value 
of  S  of  the  previous  time  step  and  compute  vt  (ut  =  c2sTt6t )  according  to  Eq.  (E.14a)  instead 
of  Eq.  (E.15).  We  will  evaluate  both  these  options.  Furthermore,  it  is  possible  to  use  finite 
difference  approximation  for  S.  We  will  also  investigate  this  option. 

It  is  important  to  point  out  the  salient  difference  between  the  LBE-LES  and  the  NS-LES. 
In  the  NS-LES,  the  eddy  viscosity  is  evaluated  and  then  used  to  determine  the  evolution 
of  the  flow  fields  in  the  next  time  step.  In  the  LBE-LES,  the  eddy  viscosity  affects  the 
relaxation  process  of  the  flow  fields  as  well  as  other  nonhydrodynamic  variables  (higher 
order  fluxes).  The  relaxation  process,  as  described  by  Eq.  (E.12),  does  not  force  the  flow 
fields  to  immediately  attain  the  expected  state  specified  by  the  equilibrium  distributions. 
This  preserves  more  spatio-temporal  memory  effects  in  the  LBE-LES  which  are  absent  in 
the  NS-LES  counterpart.  Some  preliminary  tests  using  the  LBE-LES  have  yield  encouraging 
results  (23;  25;  26;  27).  In  this  work  we  will  compare  the  LBE-LES  and  NS-LES  in  the 
fundamental  problem  of  HIT. 

In  all  the  results  presented  in  this  Appendix,  we  use  the  single-relaxation-time  LBE 
obtained  from  BGK  model  for  the  collision  operator.  Preliminary  computations  of  LBE- 
DNS  in  inertial  frame  using  Multiple-relaxation-time(MRT)  lattice  Boltzmann  model  show 
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no  distinguishable  difference  in  the  results  obtained.  We  will  investigate  the  use  of  MRT 
for  LBE-DNS  in  rotational  frame  and  LBE-LES  in  the  future  work. 


E.4  Simulation  results 


All  numerical  simulations  are  conducted  in  a  three-dimensional  periodic  cube  with  various 
resolutions:  JV3.  The  initial  incompressible  homogeneous  isotropic  velocity  field  is  generated 
in  spectral  space  with  the  following  energy  spectrum  in  a  prescribed  range  Kmjn  <  k  <  Kmax 
(For  further  details,  see  the  on-line  article:  T.  Miyauchi  and  T.  Ishizu,  “Direct  numerical 
simulation  of  HIT  decay  of  passive  scalar  fluctuation,”  available  at  http :  //cf  d .  me .  umist . 

ac.uk/ercofold/database/test48/test48.html): 


0.038/cm  exp(— 0.14  k2), 

0, 


K  £  [Kmim  Kmax]i 
K  ^  [Kmin>  ftmax]> 


(E.16) 


with  random  phase  and  then  transferred  to  physical  space.  In  the  above,  m  is  set  to  4  or  2 
for  different  simulations. 


The  initial  density  fluctuation  Sp  (or  the  pressure  p)  can  be  consistently  obtained  by 
an  iteration  procedure.  It  is  important  to  stress  that  the  preparation  of  the  initial  data  is 
crucial  in  the  LBE-DNS  simulations  of  HIT.  The  pressure  obtained  by  solving  the  Poisson 
equation  from  the  initial  velocity  u0  is  inconsistent  and  insufficient  to  initialize  the  LBE 
simulation.  It  is  inconsistent  because  LBE  is  intrinsically  compressible  thus  the  Poisson 
equation  is  not  satisfied  exactly.  It  is  insufficient  because  LBE  initial  data  consists  of  more 
than  the  hydrodynamic  variables  and  the  nonhydrodynamic  variables  cannot  be  specified 
by  solving  hydrodynamic  equations.  By  consistently  constructing  the  initial  data  for  the 
LBE  simulation,  we  are  able  to  minimize  the  error  due  to  initialization  and  the  difference 
between  the  LBE  and  pseudo-spectral  simulations  (8). 
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E.4.1  LBE-DNS  of  decaying  isotropic  turbulence 

We  first  present  the  results  from  LBE-DNS  of  decaying  HIT  in  the  inertial  frame  at  two 
resolutions:  643  and  1283.  All  initial  spectra  are  given  by  Eq.  E.16  in  which  m  is  set  to  4 
unless  indicated  otherwise. 

Figure  E.l  shows  the  evolutions  of  the  normalized  kinetic  energy  k/k$  and  the  normalized 
dissipation  rate  e/e0  with  respect  to  normalized  time  t'  =  teo/fco  for  the  cases  of  1283  and 
643.  The  parameters  for  both  cases  are  ums  =  0.023  and  v  «  0.0017  (r  =  0.505).  In  the  case 
of  643,  the  initial  energy  spectrum  is  non-zero  in  the  range  4  <  k  <  8,  resulting  in  ReA  «  53. 
For  the  case  of  1283,  the  initial  energy  spectrum  is  non-zero  in  the  range  1  <  k  <  8, 
resulting  in  ReA  ~  119.  In  the  absence  of  production,  kinetic  energy  decays  monotonically 
in  time,  whereas  at  early  stages  the  dissipation  rate  increases.  This  increase  in  e/e0  is 
completely  consistent  with  known  turbulence  physics  (explained  further  below)  and  the  same 
phenomenon  is  also  seen  in  NS-DNS  results.  Following  this  period  of  increasing  dissipation, 
both  the  kinetic  energy  and  dissipation  decay  monotonically.  The  decay  exponent  n  of  the 
kinetic  energy  in  these  low  ReA  simulations  varies  in  time.  Furthermore,  ReA  itself  is  a 
function  of  time  as  the  turbulence  decays.  The  variation  of  n  vs.  ReA  in  various  simulations 
are  shown  in  Fig.  E.2.  The  dependence  of  n  on  ReA  obtained  by  the  LBE-DNS  is  very  similar 
to  that  observed  in  NS-DNS  calculations  (3).  The  values  of  n  obtained  in  the  present  work 
agree  well  with  the  experimental  and  numerical  NS-DNS  data. 

In  Figure  E.3,  the  left  plot  shows  the  compensated  energy  spectra  [ E(k ,  i)/«4]  of  the 
above  1283  simulation  at  early  times  during  which  cascade  is  the  dominant  process.  Initially, 
the  spectrum  (dashed  line)  is  narrow  and  soon  the  energy  spreads  to  higher  wave  numbers 
(smaller  scales)  due  to  the  nonlinear  cascade  process.  This  phenomena  leads  to  the  increase 
of  the  dissipation  rate  in  physical  space,  as  shown  in  Fig.  E.l.  This  fact,  in  itself,  is 
significant  since  advection  (the  source  of  nonlinearity)  is  handled  very  differently  in  LBE. 
At  this  stage,  the  spectrum  scales  as  E(k,  t*)  ~  k4  at  small  k.  The  right  plot  shows  another 
1283  simulation  which  has  the  same  rms  of  the  initial  velocity  field  and  v  as  the  left  case 
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but  the  non-zero  energy  is  concentrated  in  the  range  of  8  <  k  <  16,  resulting  in  Rex  «  67. 
In  spite  of  the  different  range  of  initial  energy,  the  spectrum  still  scales  as  E(K,t')  ~  k4. 

Next  we  show  results  from  a  set  of  simulations  in  which  the  initial  spectrum  is  also 
given  by  Eq.  E.16  but  with  m  =  2.  In  Figure  E.4,  the  compensated  spectrum  E(k)/k2 
is  shown  at  various  times.  It  is  seen  that  the  spectrum  now  scales  as  k2.  In  summary,  as 
shown  in  Figure  E.5,  the  low  wave-number  spectra  scale  as  E(k)  ~  k4  if  m  —  4(left  plot) 
and  E(k)  ~  «2  if  m  =  2(right  plot).  This  dependence  of  low-wave  number  scaling  on  initial 
spectrum  is  in  exact  agrement  with  the  results  reported  in  (2;  3) 

Rotating  Reference  Frame.  LBE-DNS  of  decaying  HIT  in  a  rotating  frame  is  also 
performed.  Without  loss  of  generality,  we  assume  that  the  frame  of  reference  rotates  about 
the  2-axis  with  the  angular  velocity  U  =  (0,0,w).  The  Rossby  number  is  defined  as  Ro  = 
kpu rms/w>  where  kp  characterizes  the  energy  containing  wave  number  at  t  =  0.  Here,  we 
USe  Kp  =  (/Cmax  Kmin)/2- 

The  effects  of  rotation  are  scale  dependent  and  they  are  enhanced  by  increasing  the 
rotation  rate  w  (decreasing  Ro).  In  general,  it  has  been  well  understood  that  rotation  slows 
down  the  cascade  and  delays  the  approach  to  equipartition  (28;  29).  These  features  are 
captured  in  Figs.  E.6  and  E.7.  Figure  E.6  shows  the  evolution  of  kinetic  energy  at  various 
Rossby  numbers  in  a  simulation  with  1283  resolution.  The  initial  energy  spectrum  is  non¬ 
zero  in  the  range  of  1  <  «  <  8.  As  expected,  the  energy  decay  slows  down  with  decreasing 
Rossby  number  (or  increasing  rate  of  rotation).  Closer  examination  of  the  spectra  (Fig.  E.7) 
shows  the  tendency  to  maintain  more  energy  at  the  small  wave  numbers  (large  scales)  when 
the  system  rotates.  The  faster  the  system  rotates  (smaller  Rossby  number),  the  more 
prominent  is  this  tendency. 

E.4.2  LBE-LES  of  decaying  isotropic  turbulence 

In  the  previous  section,  it  was  clearly  demonstrated  that  the  LBE  method  is  an  accurate 
DNS  tool  for  turbulence.  It  is  also  important  to  asses  the  ability  of  LBE  in  the  LES  context. 


E.4.  Simulation  results 


124 


In  this  work,  we  conduct  LES  of  decaying  HIT  without  rotation.  Investigation  of  the  LES 
of  decaying  HIT  with  rotation  will  be  presented  in  the  near  future. 

In  order  to  perform  close  comparisons  with  DNS  results,  we  perform  LES  with  the  initial 
large-scales  identical  to  that  of  1283  LBE-DNS  case  (corresponding  to  the  results  presented 
in  Fig.  E.l).  Thus  the  initial  flow  fields  obtained  in  the  LBE-DNS  are  appropriately  trun¬ 
cated  in  spectral  space  to  yield  the  initial  fields  for  the  LBE-LES  for  323  and  643  resolutions. 

In  other  words,  the  initial  LBE-LES  field  is  obtained  by  filtering  out  all  wave  numbers  above 
16  for  323  and  32  for  643. 

Calibration  of  Csm.  Our  first  exercise  is  to  determine  the  appropriate  Smagorinsky 
constant  for  LBE-LES.  Figure  E.8  shows  the  energy  spectra  at  some  specific  time  instant 
with  different  Smagorinsky  constant  values  for  both  323  and  643  cases.  The  instantaneous 
LES  spectra  with  resolutions  323  and  643  are  compared  against  DNS  spectrum  at  the  same 
time.  In  general,  643  performs  better  than  323  although  at  small  «  (large  scale)  region 
both  32s  and  643  spectra  agree  well  with  the  DNS  spectrum.  The  comparison  of  the  kinetic 
energy  decay  from  the  same  runs  is  shown  in  Fig.  E.9.  FVom  both  figures,  we  find  that 
Csm  =  0.1  yields  better  results  than  the  typical  value  of  Csm  =  0.17  used  in  the  NS-LES(9). 

The  need  for  a  reduced  Csm  in  LBE-LES  compared  with  the  NS-LES  value  can  be  explained 
as  follows. 

In  the  lattice  Boltzmann  equation  there  are  nonhydrodynamic  variables  which  are  higher 
order  fluxes  (Sfuy,  n  >  1).  These  higher  order  fluxes  could  lead  to  a  higher  effective  viscosity 
and  other  nonlinear  effects.  Thus  a  smaller  Csm  value  in  LBE-LES  is  adequate  to  achieve 
the  same  effect  as  a  larger  Csm  in  NS-LES.  In  Fig.  E.10,  we  compare  the  instantaneous 
flow  structure  of  uz(i,  j,  k  =  N/2,  t')  obtained  by  the  LBE-LES  with  that  by  LBE-DNS. 
As  shown  in  Fig.  E.10,  the  LBE-LES  appears  to  capture  the  flow-field  structure  quite 

adequately  even  with  a  coarse  resolution  of  323.  In  all  subsequent  calculations,  we  use 

CSm  —  0.1. 

Other  methods  of  S  computation.  In  all  the  above  calculations,  we  determine 
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S{m  the  express  of  i/t)  from  the  second  moment  of  nonequilibrium  distribution  functions  at 
the  current  time-step.  We  now  investigate  the  other  options  for  determining  S  which  may 
provide  some  computational  advantages. 

First,  we  test  two  different  ways  to  compute  the  strain  rate  tensor:  (i)from  the  second 
moment  of  nonequilibruim  distribution  functions  as  given  in  Eq.  E.ll  and  (ii)by  finite- 
difference  approximation  of  derivatives  Si:,-  =  (c^  +  ^)/2.  In  Figure  E.ll,  we  compare 
the  kinetic  energy  evolution  from  both  computations  with  LBE-DNS  result  at  two  reso¬ 
lutions.  In  the  323  case(left  plot),  it  is  seen  that  the  computation  using  finite-difference 
approximation  quickly  diverges  while  computation  from  nonequilibrium  distribution  func¬ 
tion  moment  captures  the  DNS  result.  In  the  643  case(right  plot),  since  the  resolution  is 
high  enough,  both  computations  yield  good  results.  Next,  we  test  two  different  ways  of 
computing  strain  rate  S  from  the  second  moment  of  nonequilibruim  distribution  functions. 
In  the  implicit  method,  S  of  the  current  time-step  is  used  to  yield  Eq.  E.15  for  vt.  In  the 
explicit  approach,  S  from  previous  time-step  is  used  (as  in  Eq.  E.14a).  Figure  E.12  depicts 
the  contours  of  the  instantaneous  flow  structure  of  uz{i,  j,  k  =  N/ 2,  t')  obtained  from  LBE- 
LES  (323)  by  using  the  two  formulae  to  compute  vt.  The  velocity  fields  obtained  with  the 
two  formulae  are  almost  identical,  as  shown  in  Fig.  E.12;  the  L2-norm  difference  between 
the  two  velocity  fields  is  less  than  0.02%.  Therefore,  we  verify  that  the  eddy  viscosity  can 
be  computed  by  either  Eq.  (E.15)  or  Eq.  (E.14a)  without  significant  effect  on  the  flow  fields. 

E.4.3  LBE-LES  vs.  NS-LES 

We  further  compare  the  LBE-LES  with  the  NS-LES  results  at  323  resolution  and  the  re¬ 
sults  are  shown  m  Figs.  E.13  -  E.15.  The  initial  velocity  fields  for  LBE-LES  and  NS-LES 
calculation  are  nearly  identical.  FVom  Fig.  E.13  and  E.14,  it  is  seen  that  the  kinetic  energy 
and  spectra  computed  from  LBE-LES  are  somewhat  closer  to  the  DNS  results  than  those 
calculated  from  NS-LES.  The  difference  can  be  seen  more  clearly  from  contours  of  the  in¬ 
stantaneous  flow  field  uz(i,  j,  k  =  N/ 2,  t')  at  three  different  times  as  shown  in  Fig.  E.15. 
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By  Comparing  the  LBE-LES  results  (left  column)  and  the  NS-LES  results  (right  column) 
with  the  LBE-DNS  results  (center  column),  we  observe  that  LBE-LES  preserves  the  flow 
structure  better  than  the  corresponding  NS-LES.  The  323  LBE-DNS  contours  shown  here 
are  obtained  by  truncating  the  1283  LBE-DNS  data. 

E.5  Summary  and  conclusions 

In  this  Appendix  we  perform  DNS  and  LES  of  the  classical  decaying  HIT  problem  with 
and  without  reference  frame  rotation  using  the  LBM.  Three  categories  of  simulations  have 
been  performed.  First  we  conduct  the  direct  numerical  simulations  by  using  the  LBM. 
Well  known  power-law  decay  of  the  kinetic  energy  is  reproduced.  The  decay  exponents 
obtained  in  the  LBE  simulations  are  in  good  agreement  with  the  results  from  experimental 
measurements  and  NS-DNS  calculations.  The  low-wavenumber  energy  spectrum  scaling 
depends  on  initial  conditions.  Both  and  k2  scaling  are  obtained  from  appropriate  initial 
conditions  consistent  with  (2;  3).  The  effect  of  rotation  on  turbulence,  which  is  to  suppress 
the  spectral  cascade,  is  also  well  captured. 

Second,  we  conduct  a  comparative  study  of  the  LBE-LES  and  the  LBE-DNS.  Compar¬ 
isons  between  643,  323  LBE-LES  and  1283  LBE-DNS  show  that  the  large  scale  motion  is 
well  captured  by  LBE-LES.  A  smaller  Smagorinsky  constant,  Csm  =  0.1,  is  demonstrated 
to  achieve  better  performance  in  LBE-LES.  By  choosing  appropriate  Smagorinsky  con¬ 
stant,  even  323  can  adequately  capture  large  scale  motions.  Finally  we  compare  both  the 
LBE-LES  and  NS-LES  with  the  corresponding  DNS  results.  Our  comparisons  indicate  that 
LBE-LES  has  better  capability  to  preserve  flow-field  structure  than  NS-LES.  This  work 
further  establishes  the  LBM  as  a  viable  computational  tool  for  turbulence  simulations. 
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Figure  E.l:  Time  evolution  of  the  normalized  kinetic  energy  k/k0  (solid  lines)  and  normal¬ 
ized  dissipation  rate  e/eo  (dashed  lines)  for  643  (lines  and  symbols)  and  1283  (lines  only) 
by  using  LBE-DNS. 
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Figure  E.2.  Dependence  of  the  decay  exponent  n  —  \/{C2  —  1)  on  initial  conditions  and 
Re*.  The  quantity  C2  is  depicted  in  the  figure  in  stead  of  n.  Solid  lines  represent  NS-DNS 
data  from  (3)  and  symbols  correspond  to  the  LBE-DNS  results  of  the  present  work.  For 
the  1283  resolution,  •:  urms  =  0.0064,  kmili  =  1,  kmax  =  8,  and  v  =  0.01  (r  =  0.53);  A: 
urms  =  0.021,  kmin  =  8,  kmax  =  16,  and  v  «  0.00167  (r  =  0.505);  o:  nrms  =  0.022,  kmin  =  1, 
kmax  =  8,  and  v  «  0.00167  (r  =  0.505).  For  the  643  resolution  (x):  urms  =  0.022,  fcmin  =  4 
and  kmax  =  8,  and  v  «  0.00167  (r  =  0.505). 
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Figure  E.3:  Compensated  energy  spectra  for  two  cases  of  1283  at  early  times,  t'  =  0.022, 
0.044,  0.066,  and  0.088(left)  and  t'  =  0.0022,  0.022,  and  0.088(right).  The  dashed  lines 
represent  the  initial  spectra  given  by  Eq.  (E.16)  with  m  =  4.  The  spectra  scale  as  E(k,  t')  ~ 
k4  at  small  k 
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Figure  E.4:  Compensated  energy  spectra  for  a  1283  simulation  (urms  =  0.023,  u  a  0.0017 
(r  =  0.505),  and  Re\  a  141),  at  early  times,  t'  =  0.011,  0.017,  0.027.  The  dashed  lines 
represent  the  initial  spectra  given  by  Eq.  (E.16)  with  m  =  2.  The  spectra  scale  as  E(k,  t')  ~ 
k2  at  small  k. 
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Figure  E.5:  Energy  spectra  for  two  cases  of  1283  above  with  initial  energy  concentrat¬ 
ing  in  the  range  of  1  <  k  <  8  at  t'  —  0  and  t'  =  0.022  respectively.  The  ini¬ 
tial  energy  spectra  are  given  by  E(n)  =  0.038k4  exp(— 0.14  k2)  (upper  plot)  and  E(k)  = 
0.038k2  exp(— 0.14  K2)(bottom  plot)  respectively. 
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Figure  E.6:  Kinetic  energy  decay  in  1283  LBE-DNS  with  different  Rossby  number  Ro. 


Figure  E.7:  Energy  spectra  at  t'  =  10.5  with  different  Rossby  number  Ro  for  643  case: 
urms  =  0.023,  fcmir.  =  1,  fcmax  =  4,  and  v  =  0.01  (r  =  0.53).  The  dashed  line  is  the  inertial 
case  (fl  =  0  or  Ro  =  oo). 
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Figure  E.8:  Energy  spectra  at  t'  =  0.04796  with  different  Csm  and  resolution.  Solid  line  for 
LBE-DNS  (1283)  and  symbols  for  LBE-LES  (o:  643,  Csm  =  0.1;  A:  643,  Csm  =  0.17;  x: 
323,  Csm  =  0.1;  +:  323,  Csm  =  0.17). 
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Figure  E.9:  Kinetic  energy  decay  with  different  Carn  and  resolution.  Solid  line  for  LBE- 
DNS  (1283)  and  symbols  for  LBE-LES  (o:  643,  Csm  =  0.1;  A:  643,  Csm  =  0.17;  x:  32s, 
Csm  =  0.1;  +:  323,  Csm  =  0.17). 
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DNS,  323  DNS,  643 


Figure  E.10:  Contours  of  the  instantaneous  flow  field  uz(i,  j,  k  =  N/2,  t').  LBE-DNS  vs. 
LBE-LES  with  different  resolutions.  The  323and  643  LBE-DNS  contours  shown  here  are 
obtained  by  truncating  the  1283  LBE-DNS  data 
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Figure  E.12:  Contours  of  the  instantaneous  flow  field  uz(i,  j,  k  =  N/2,  t')  obtained  by 
LBE-LES  with  a  resolution  of  323  and  two  different  formulae  for  vt  in  two  different  times. 
Top  row:  ft,  is  computed  according  to  Eq.  (E.15)  and  bottom  row:  is  computed  according 
to  Eq.  (E.14a)  with  one  time  step  lagging  in  S. 
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Figure  E.13:  Comparison  of  kinetic  energy  decay  of  323  LBE-LES(o),  323  NS-LES(*)  and 
1283  LBE-DNS(A). 
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LBM-LES,  t'=0.006  LBM-DNS,  t'=0.006  NS-LES,  t'=0.006 


LBM-LES,  t’ =0.01 5  LBM-DNS,  t'=0.01 5  NS-LES,  t'=0.015 


Figure  E.15:  Contours  of  the  instantaneous  flow  field  uz(i,  j,  k  =  N/2,  t').  The  LBE-LES 
and  NS-LES  with  a  resolution  of  323  compared  to  the  LBE-DNS  with  a  resolution  of  1283 
at  three  different  times.  The  323  LBE-DNS  contours  shown  here  are  obtained  by  truncating 
the  1283  LBE-DNS  data 
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Appendix  F 


Study  of  Homogenous  Turbulence 
Shear 


Nomenclature 

bij 

Anisotropy  tensor 

<  UiUj  > 

Reynolds  stress 

Pij 

Production  tensor 

Rij 

Pressure-strain  tensor 

e 

Dissipation  rate 

k 

Turbulent  kinetic  energy 

Dissipation  tensor 

Rex 

Taylor  microscope  Reynolds  num¬ 
ber 

p 

Production  of  k 

UJ 

Frequency  of  shear 

s 

Shear  rate 

Sk/e 

Normalized  shear  rate 

cs 

Speed  of  sound 

ea 

Discretized  particle  velocity 

f 

Distribution  function 

Equilibrium  distribution  function 

U 

Fluid  velocity 

P 

Fluid  density 

St 

Time  step 

sx 

Grid  size 

T 

V 

Normalized  relaxation  time 
Viscosity 

wa 

Weight  function 

F.l  Introduction 


In  Appendix  E,  we  have  studied  decay  isotropic  turbulence  in  inertial  and  rotational  frame. 
As  a  part  of  our  effort  to  perform  LBM-DNS  of  turbulence,  we  will  investigate  the  homoge- 
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nous  shear  turbulence  in  this  Appendix. 

Homogenous  shear  flow  has  been  extensively  studied  by  experiments  and  numerical 
simulations.  Rose  (1),  Champagne  Harris  and  Corrsin  (2),  Tavoularis  and  Corrsin  (3;  4) 
studied  homogenous  shear  turbulence  in  tunnel  to  investigate  the  local  isotropy  of  small 
scales.  The  issue  of  local  isotropy  was  also  addressed  numerically  by  Pumir  and  Shraiman 
(5)  and  Pumir  (6).  Rogers  and  Moin  examined  the  structure  of  the  vorticity  in  homogenous 
shear  flows  (7).  They  confirmed  the  presence  of  hairpin  vortices  in  homogeneous  shear  flow. 
Lee  et  al  (8)  compared  the  structure  of  turbulence  under  high  shear  rate  in  shear  flow  with 
that  in  turbulent  channel  flow  and  found  many  dynamical  features  in  the  two  flows  to  be 
similar.  Kida  and  Tanaka  (9)  discussed  the  regeneration  cycle  of  the  streamwise  vortices 
in  a  homogeneous  shear  flow.  Jacobitz  et  al.  (10)  studied  uniformly  sheared  and  stably 
stratified  flows.  They  found  that  the  turbulence  evolution  depends  strongly  on  the  initial 
value  of  the  Taylor  microscale  Reynolds  number  and  the  initial  value  of  the  normalized 
shear  number  Sk/e  in  both  cases.  Most  of  the  above  simulations  used  the  pseudo-spectral 
method  of  the  solution  of  Navier-Stokes  equations  (11). 

In  this  Appendix,  we  will  first  verify  the  reliability  of  LBM  in  the  DNS  of  homogenous 
turbulence  and  then  focus  on  the  new  physics  in  homogenous  turbulence  flow  subjected 
to  periodically  varied  shear.  This  Appendix  is  structured  as  follows.  The  Reynolds  stress 
evolution  equations  are  given  in  Section  F.2.  In  Section  F.3,  firstly,  the  classical  results 
(3;  4;  7)  are  used  to  verify  the  viability  of  LBM  in  DNS  of  homogenous  turbulent  shear 
(constant  shear).  Secondly,  results  of  LBM-DNS  of  homogenous  turbulence  subjected  to 
periodically  shear  are  presented  and  the  new  physics  in  this  flow  are  discussed.  Finally, 
we  compare  the  DNS  results  with  the  results  obtained  by  Reynolds  average  Navier-stokes 
(RANS)  turbulence  models  in  the  periodically  varied  shear  cases.  These  analysis  provide 
some  new  knowledge  for  developing  models  in  turbulence  subject  to  unsteady  forcing.  A 
brief  conclusion  is  presented  in  Section  F.4. 
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F.2  Reynolds-stress  evolution  equations 

The  Reynolds-stress  evolution  equation  can  be  written  as 

D  .  a 

Jyj.  \uiuj)  "I"  —  Pij  +  Rij  ~~ 

where  the  turbulent  transport  of  Reynolds-stress  ( T^j )  can  be  decomposed  as 

7\.  .  —  rpiu)  i  rp(p)  ,  rp(v) 

±ki3  “  1kij  +  J-kij  +  1kij- 


(P.1) 


In  the  above, 


Tkij  -  (uiujuk)  ,  =  ~  (UiP  )  ^  (^jp) 

Further,  production,  pressure  redistribution,  and  dissipation  are  respectively: 

D  _  /  x  d(Uj)  ,  v  d(Ui ) 

pij  =  -  {um)  ~  (ujuk)  -j-L, 

/p  f  dui  duj 
j  \  p  V&r,  dxi 


'«■>(££) 


(F.2) 

(F.3) 

(F.4) 


where  ()  represents  the  mean  of  the  fluctuating  variables;  U  is  fluid  velocity;  u  and  p  are 
fluctuating  velocity  and  pressure;  and  p  is  the  fluid  density. 

In  homogenous  turbulence,  the  evolution  equation  simplifies  to 


—  (uzUj)  =  +  Rij  —  fyj. 


(F.5) 


In  our  simulation,  the  shear  is  only  applied  in  the  flow  along  y-direction.  For  pertinent 
components  of  Reynolds  stress  in  this  homogenous  turbulence  shear  flow,  we  get 


d_ 

at 


d(Ui)  /Jdui 

p  dx\ 


_ <„,„,)  .  -  (U1„2)  (P.6) 


dui  dui  \ 


9  ,  .  i  p  du2  \  /  du2  du2 


(F.7) 
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~  <W3)  =  - 


yP_  dU'j 

"  p  dx3 


-2v 


du3  dus  \ 

dx  k  dxk  J 


(F.8) 
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In  studying  homogenous  turbulence  shear,  it  is  more  common  to  present  results  in  the 
anisotropy  tenser  which  is  defined  as  bij  =  —  ^5{j. 


F.3  Results  and  discussions 

F.3.1  Homogenous  shear  flow 

As  mentioned  in  introduction,  homogenous  turbulent  shear  has  been  extensively  studied. 
In  this  section,  we  will  use  this  flow  as  a  benchmark  to  verify  the  applicability  of  LBM  in 
DNS  of  turbulence.  The  computational  domain  is  shown  in  Fig.  F.l.  The  mesh  is  uniform 
and  has  the  resolution  of  128^.  In  x-  and  z-directions,  periodical  boundary  condition  is 
used;  Non-slip  condition  is  used  at  upper  and  lower  walls.  Usually,  bounce  back  scheme  is 
used  at  the  wall.  However,  to  avoid  the  instability  in  computation,  we  use  the  equilibrium 
distribution  function  in  the  present  simulations.  For  more  details  of  discussion  of  boundary 
condition  in  LBM,  please  see  reference  (13).  In  the  simulation,  19-bit  (D3Q19)  model  is 
used  because  of  its  advantages  in  matters  of  stability  and  efficiency.  Initial  velocity  filed 
is  given  by  isotropic  turbulence.  Uniform  mean  velocity  gradient  S=||  is  imposed  on  the 
isotropic  turbulence  by  adding  linear  mean  velocity  profile  along  y-direction  at  initial  stage. 
To  maintain  the  linear  mean  velocity  profile,  the  lower  and  upper  walls  move  with  a  constant 
speed  of  U  in  opposite  direction  along  a;- axis.  At  initial  stage,  normalized  shear  Sk/e  is  3.3, 
and  Taylor  microscale  Reynolds  number  Rex  is  33.  Re\  =  2k\/v  represents  the  Reynolds 
number  based  on  the  Taylor  microscale.  Statistical  quantities  of  turbulence  are  obtained  by 
spatial  averaging  over  the  inner  128  x  64  x  128  region  in  computational  domain  to  avoid  the 
wall  effects.  The  most  important  inference  has  been  made  from  the  studies  of  Tavoularis 
and  Corrsion  (3;  4)  and  Rogers  and  Moin  (7)  is  that  the  homogenous  turbulence  reaches 


F.3.  Results  and  discussions 


148 


an  asymptotic  stage  after  initial  transition.  At  this  asymptotic  stage,  when  normalized  by 
the  imposed  shear  rate  S  and  the  kinetic  energy  k(t),  the  statistics  become  independent  of 
time  (12). 

Figure  F.2  shows  the  development  of  anisotropy  tensor  obtained  from  the  present  cal¬ 
culations  as  well  as  previous  investigations.  It  is  difficult  to  compare  present  results  with 
those  obtained  by  NS  solver  at  early  stage  because  of  difference  in  initial  conditions.  At 
the  asymptotic  stage,  it  can  be  seen  that  LBM  recovers  the  experimental  and  numerical 
results  very  well.  Fig.  F.3  compares  present  &12  with  the  results  obtained  by  Jacobitz  et  al. 
(10).  As  is  clearly  seen,  the  agreement  is  again  good.  Fig.  F.4  shows  the  development  of 
normalized  shear  and  the  production  over  dissipation  ratio  {P/e).  We  see  that  after  initial 
period,  normalized  shear  and  P/e  almost  maintain  constant  values  which  are  3.8  for  nor¬ 
malized  shear  and  1.2  for  P/e.  As  shown  by  Jacobitz  et  al  (10),  these  values  are  depended 
on  the  initial  shear  and  Re\. 

From  the  above  results,  we  can  conclude  that  the  LBM  is  an  effective  tool  for  simulating 
turbulence. 

F .3.2  Homogenous  turbulence  subjected  to  periodically  varied  shears. 

In  the  previous  case,  shear  was  constant  which  is  the  classical  problem.  Here  we  let  shear 
vary  with  time  and  investigate  the  effects  of  the  variation  on  the  development  of  turbulence. 
In  simulations,  all  the  conditions  are  kept  the  same  as  those  in  the  constant  shear  case  except 
that  non-penetration  condition  is  used  at  upper  and  low  boundaries.  Further,  shear  is  varied 
with  time  according  to  sin(wt),  where  §|  equals  to  the  one  in  constant  shear  case.  The 
variation  of  shear  is  enforced  by  adding  body  force  in  the  lattice  Boltzmann  Equation  (14), 

fa{xi  +  eaSt,t  +  6t)  -  fa(xi,t)  =  -~[fa{xi,t)  -  f^q\xht)]  -+-  5tFa 

T 

where  Fa  —  (1  —  ^p)wa  •  F,  F  =  U (y)ui  cos(u)t)  •  i,  and  i  is  the  unit  vector 

along  x-direction.  Frequency  of  the  variation  of  shear  (a;)  is  determined  based  on  the 
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magnitude  of  S.  In  simulations,  two  different  ratios  of  a;  to  5  are  used  which  are  ^=1  and 
0.5,  respectively. 

Figures  F.5  and  F.6  show  the  evolution  of  turbulent  kinetic  energy  for  the  low  frequency 
case  (|=0.5)  and  high  frequency  case  (g=1.0),  respectively.  Three  important  observations 
are  made  in  these  flows.  First,  there  exists  a  critical  frequency  at  which  turbulence  changes 
from  growth  model  to  decay  mode.  As  can  be  clearly  seen  that,  in  low  frequency  case, 
k  increases,  while  in  high  frequency  case,  k  decays.  The  second  phenomenon  is  that  the 
frequency  of  the  variation  of  k  is  twice  as  that  of  shear  in  both  cases.  The  third  phenomenon 
is  that  in  high  frequency  case,  the  variation  of  k  is  not  as  that  in  the  low  frequency  case. 
In  low  frequency  case,  the  variation  of  k  can  be  described  by  a  single  frequency  function. 
In  high  frequency  case,  k  has  two  peaks  with  different  magnitudes  during  a  period  of  shear. 
These  phenomena  will  be  explained  in  following  discussions  by  analyzing  the  variation  of 
shear  and  612  because  these  two  terms  determine  the  production  of  turbulence. 

We  will  first  look  at  the  low  frequency  case.  The  variation  of  shear  and  612  are  shown  in 
Fig.  F.7.  Just  before  time  A,  the  shear  is  negative  and  612  is  positive.  Thus  production  for 
k  is  positive.  After  time  A,  shear  becomes  positive  while  612  still  maintains  positive  from 
time  A  to  B.  Thus,  negative  production  is  generated  for  k  from  time  A  to  B  as  shown  in  Fig. 
F.8  and  is  reflected  also  in  the  rapid  decay  of  k  in  the  same  period  of  time  in  Fig.  F.5.  After 
shear  taking  effect,  612  rapidly  decreases  and  becomes  negative.  Then  production  turns  to 
be  positive  from  time  B  to  C  in  Fig.  F.8  and  is  reflected  also  in  the  growth  of  k  in  same 
time  in  Fig.  F.5.  In  the  constant  shear  case,  after  reaching  its  largest  magnitude,  612  will 
almost  be  maintained  at  that  value.  In  Fig.  F.7,  we  can  see  that  612  has  possibly  reached 
the  maximum  value  because  the  magnitude  of  612  has  reached  around  0.2  (compared  with 
0.18  in  constant  shear  case)  and  it  does  not  vary  too  much  for  a  long  time.  We  will  refer 
that  maximum  value  of  £>12  as  saturation  value. 

In  Fig.  F.8,  from  time  A  to  B,  the  turbulence  is  adjusting  to  mean  flow  and  P/e  is 
negative.  The  largest  negative  P/e  happens  somewhere  between  time  A  and  B.  It  can  be 
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seen  in  Fig.  F.7  that  shear  and  612  are  relative  small  in  period  AB.  So  the  magnitude  of 
negative  P/e  is  small.  From  time  B  to  C,  Turbulence  has  adjusted  to  mean  flow  and  P/e 
is  positive.  Shown  in  Fig.  F.7,  Concurrence  of  large  shear  and  bn  produces  larger  positive 
P/e  during  time  BC.  The  ratio  of  period  AB  to  BC  is  about  2:3,  so  turbulence  stays  in 
positive  production  much  longer  than  it  is  in  negative  production. 

Figures  F.9  shows  the  variation  of  612  when  frequency  is  high.  The  difference  between 
high  and  low  frequency  cases  is  obvious.  The  first  is  the  asymmetry  of  612  to  x-axis  in  the 
high  frequency  case:  the  magnitude  of  positive  612  is  much  smaller  than  that  of  negative 
612.  This  is  caused  by  the  initial  moving  direction  of  shear  and  the  difference  of  frequency. 
In  our  simulation,  the  shear  turns  to  positive  first.  Thus  b\ 2  becomes  negative  first  and 
may  reach  largest  negative  value.  When  shear  becomes  negative,  612  needs  to  be  turned 
to  positive.  However,  there  is  no  enough  time  for  612  to  fully  attain  the  largest  positive 
value.  Thus,  the  asymmetry  of  b\ 2  is  generated.  As  shown  in  Fig.  F.9,  the  asymmetry  of 
612  diminishes  with  time.  The  second  difference  between  high  frequency  and  low  frequency 
cases  is  the  relative  position  of  shear  and  612.  Examining  in  one  period  in  Fig.  F.7  and 
F.9b,  it  can  be  seen  that,  from  time  A  to  C,  the  relative  positions  of  shear  and  b\ 2  are 
same  for  both  cases.  This  relation  gives  a  small  negative  production  first  then  follows  a 
large  positive  production.  However  from  time  C  to  E,  the  difference  appears.  In  the  high 
frequency  case  period  CD  is  much  longer  and  period  DE  is  much  shorter  than  they  are  in 
low  frequency  case.  This  relation  generates  a  large  negative  production  during  time  CD  and 
small  positive  production  during  time  DE  in  high  frequency  case.  This  also  explains  the 
two  different  peaks  in  the  variation  of  k  in  Fig.  F.6a.  For  high  frequency  case,  the  negative 
and  positive  production  will  offset  each  other  and  produce  no  net  effect  in  the  development 
of  turbulence. 

It  can  be  seen  that  magnitude  of  P/s  increases  with  time  in  Fig.  F.10  for  high  frequency 
case.  This  is  primarily  due  to  the  decrease  of  e.  The  Mean  value  of  P/e  is  just  above  zero 
and  much  less  than  1.  From  fig.  F.10,  we  can  see  that  the  decay  of  turbulence  in  high 
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frequency  is  not  caused  by  the  small  magnitude  of  P/e.  It  is  caused  by  large  negative  P/e. 
Negative  production  means  turbulence  transfers  its  kinetic  energy  into  mean  flow.  The 
mechanism  of  this  phenomenon  will  be  investigated  in  details  in  future  works. 

The  details  of  development  of  turbulence  anisotropy  tensor  are  given  in  Figs.  F.ll  and 
F.12.  In  low  frequency  case,  anisotropy  in  turbulence  has  been  well  developed.  The  mean 
values  of  and  622  3-re  about  0.15  and  -0.1  respectively.  In  high  frequency  case,  mean 
values  of  anisotropy  &11  and  622  are  about  0.  Thus,  turbulence  has  no  opportunity  to  develop 
in  high  frequency  case. 

To  gain  more  insight  of  evolution  of  turbulence,  We  investigate  the  terms  in  the  Reynolds- 
stress  equation  (Eqs.  F.6-F.9).  Figs.F.13-F.14  show  the  comparison  of  time  derivative  of 
Reynolds-stress.  We  use  two  ways  to  obtain  those  values.  One  is  from  the  right  hand  side  of 
Reynolds-stress  evolution  equation  (Eqs.  F.6-F.9)  and  the  other  is  from  the  finite  difference 
operation  on  Reynolds-stress: 

d  (uiUj)  _  ( UiUj(t  +  dt ))  -  ( UiUj(t  -  dt)) 

dt  ~  26t 

The  results  from  two  different  methods  agree  very  well.  This  guarantees  that  each  term  in 
the  Reynolds-stress  equation  budget  is  reliable. 

Figure.F.15  shows  the  different  terms  in  Reynolds  stress  equation  for  {u\u2)  (Eq.  F.9). 
The  pressure  redistribution  and  dissipation  always  counter  the  effect  of  production.  The 
magnitude  of  dissipation  is  smaller  than  those  of  production  and  redistribution.  Fig.  18 
shows  the  budget  of  (uiui).  Pressure  redistribution  almost  always  counters  the  effect  of 
production  except  when  production  is  close  to  zero.  In  the  earlier  stage,  magnitude  of 
dissipation  is  comparable  to  those  of  production  and  redistribution.  In  the  later  stage,  it  is 
very  small.  In  the  Reynolds  equations  for  (u2u2)  and  (u3u3),  the  production  is  zero.  The 
periodic  variation  of  dissipation  is  relatively  small.  Pressure  redistribution  plays  a  key  role 
in  determining  the  variation  of  kinetic  energy  as  shown  in  Figs.  F.17  and  F.18.  In  Fig. 
F.19,  the  variations  of  pressure  -strain  are  shown.  It  can  be  seen  that  pressure  redistribution 
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terms  Rll  and  R22  evolve  much  similarly,  while  phase  of  R33  is  almost  180  degree  different 
from  those  of  Rll  and  R22.  This  means  that  both  Rll  and  R22  will  transfer  turbulence 
kinetic  energy  from  (uiui)  and  (U2U2)  to  (11,3113).  This  is  different  from  wall  turbulence  flow. 
In  the  wall  flow,  R22  will  help  to  maintain  turbulence  energy  in  (u2U2)-  As  is  well  known, 
the  net  effect  of  pressure  redistribution  is  zero  which  is  verified  in  Fig.  F.19. 

F.3.3  Comparison:  DNS  vs.  RANS  models 

In  this  Section,  we  compare  DNS  results  with  those  obtained  by  solving  Reynolds  averaged 
Navier-Stokes  equations  (RANS)  with  turbulence  models.  Here,  we  use  three  different 
models  which  are  LRR-IP,  LRR-QI,  and  SSG  models  (12).  In  the  high  frequency  case,  DNS 
and  models  predict  decay  of  A;  as  shown  in  Fig.  F.20.  But  models  predict  slow  decay  rate.  In 
low  frequency  case,  models  fail  to  predict  growth  of  k  as  shown  in  Fig.  F.21.  Comparisons 
of  results  of  &12  are  given  in  Figs.  F.22  and  F.23.  The  agreement  is  somewhat  better.  This 
is  to  be  expected  as  &12  evolution  is  dominated  by  shear.  However,  in  both  cases,  model 
yields  smaller  variations  of  612  than  DNS.  In  the  low  frequency  case,  612  does  not  reach  the 
same  saturation  values  as  in  DNS.  In  high  frequency  case,  models  predict  very  large  values 
of  6n.  In  low  frequency  case,  SSG  model  gives  the  closest  result  to  that  of  DNS  for  bn- 

F.4  Conclusions 

DNS  of  homogenous  turbulence  subjected  to  constant  shear  is  performed  using  Lattice 
Boltzmann  method.  The  results  agree  very  well  with  classical  results  obtained  by  using 
Navier-Stokes  equations.  Suitability  of  lattice  Boltmann  method  for  DNS  of  turbulence  has 
been  thus  verified.  In  DNS  of  homogenous  turbulence  subjected  to  periodically  varied  shear, 
a  critical  frequency,  which  determines  whether  k  decays  or  not,  is  found.  It  seems  that  612 
reaching  a  saturation  value  is  a  necessary  condition  in  order  to  maintain  turbulence.  We 
observe  negative  production  for  k  which  means  there  is  a  mechanism  to  take  kinetic  energy 
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away  from  turbulence  and  transfer  it  to  mean  flow.  The  details  of  this  mechanism  need  to 
be  studied  further.  The  comparisons  of  results  of  DNS  and  turbulence  models  show  that 
models  fail  to  predict  the  critical  frequency.  In  terms  of  anisotropy  of  turbulence,  models 
perform  a  little  better. 


r 
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Figure  F.2:  The  results  of  anisotropic  tensor  obtained  from  DNS  using  LBM,  Navier-Stokes 
equations. 


Figure  F.3:  The  results  of  bn  obtained  from  DNS  using  LBM  and  Navier-Stokes  Equations 
by  Jacobitz  et  al.  (10)  with  the  initial  values  of  ReA  =  44.72  and  SK/e  =  2.0. 
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Figure  F .4:  Dimensionless  shear  rate  and  production  over  dissipation  ratio  by  DNS  using 
LBM.  Solid  circle  is  the  value  of  shear  obtained  in  the  DNS  by  Rogers  and  Moin  (7). 


Figure  F.5:  Evolution  of  k  for  ^  —  0.5  case 
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F.5.  Figures 


d<u22u22>/dt  using  numerical  diffrential 
d<u22u22>/dt  using  Reynold  stress  Eqn 

<u22u22>  is  normalized  by  <u22u22>  att=0 


Figure  F.14:  Comparison  of  time  derivative  of  (1^2)  for  %  =  1.0  case 
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Figure  F.15:  Variations  of  terms  in  Reynolds-stress  equation  for  (uiu2)  for  f  =  1.0  case 
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(a) 


(b) 

Figure  F.16:  (a)Variations  of  terms  in  Reynolds-stress  equation  for  (ui«i)  at  earlier  stage 
for  |  =  l.Ocase;  (b)  Variations  of  terms  in  Reynolds-stress  Equation  for  (uiui)  at  later 
stage  for  |  =  1.0  case 
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Figure  F.17:  Variations  of  terms  in  Reynolds-stress  equation  for  (u2u2)  at  later  stage  for 
5  =  1.0  case 


Figure  F.18: 
=  1.0  case 


Variations  of  terms  in  Reynolds-stress  equation  for  («3u3)  at  later  stage  for 
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Figure  F.21:  Comparison  of  k  obtained  by  DNS  and  Reynolds  averaged  Navier-Stokes 
turbulence  models  for  ^  =  0.5  case 


DNS 

LRR-IP  model 


Figure  F.22:  Comparison  of  bi 2  obtained  by  DNS  and  Reynolds  averaged  Navier-Stokes 
turbulence  models  for  ^  =  1.0  case 
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Figure  F.23:  Comparison  of  612  obtained  by  DNS  and  Reynolds  averaged  Navier-Stokes 
turbulence  models  for  fj?  =  0.5  case 


Figure  F.24:  Comparison  of  bn  obtained  by  DNS  and  Reynolds  averaged  Navier-Stokes 
turbulence  models  for  ^  =  1.0  case 
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Figure  F.25:  Comparison  of  bn  obtained  by  DNS  and  Reynolds  averaged  Navier-Stokes 
turbulence  models  for  ^  =  0.5  case 
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Appendix  G 


Lattice  Boltzmann  Models  for 
Axisymmetric  Flows  in  Cylindrical 
Coordinate 


Nomenclature 


Cs 

Speed  of  sound 

Bqc 

f 

Distribution  function 

f(eq) 

u 

Fluid  velocity 

p 

St 

Time  step 

Sx 

T 

Normalized  relaxation  time 

wa 

Po 

Mean  density 

P' 

E2 

L2  norm  error 

h 

V 

Viscosity 

Discretized  particle  velocity 

Equilibrium  distribution  function 

Fluid  density 

Grid  size 

Weight  function 

Density  fluctuation 

Grid  size 


G.l  Introduction 

LBM  has  demonstrated  a  significant  potential  and  broad  applicability  with  numerous  ad¬ 
vantages  both  physically  and  computationally,  including  the  parallelism  of  algorithm,  the 
simplicity  of  programming  and  the  ability  to  incorporate  microscopic  interactions.  Several 


171 


G.2.  Lattice  Boltzmann  model  for  axisymmetric  flows 


172 


shortcomings,  however,  remained  to  be  addressed.  One  of  them,  which  is  focused  on  partly 
in  this  Appendix,  is  the  extension  of  LBM  from  the  Cartesian  coordinates  to  curvilinear 
coordinates  to  meet  the  needs  of  many  practical  flow  problems.  Several  efforts  have  been 
made  in  this  issue  and  achieved  some  successes.  The  finite  volume  version  of  LBM  by 
Succi  and  his  coworkers(l;  2)  and  the  interpolation-supplemented  lattice  Boltzmann  equa¬ 
tion  model  by  He  et  al(3;  4)  extended  LBM  to  curvilinear  coordinate  system.  Yu  and 
Zhao(5)  applied  LBM  in  paraboloidal  coordinate  system  to  simulate  Rossby  vortex.  They 
specially  introduced  a  scaling  factor  to  re-scale  all  related  physical  quantities  and  to  map 
non-uniform  curvilinear  mesh  onto  uniform  flat  mesh.  All  above  methods  focus  on  applying 
current  lattice  Boltzmann  model  in  irregular  mesh.  Recently,  Halliday  et  al(6)  presented  a 
model  in  cylindrical  coordinate  for  pseudo-Cartesian  flow.  However,  the  implementation  of 
this  model  is  very  complicated  and  not  practical  somehow. 

In  this  Appendix,  we  present  a  lattice  Boltzmann  model  in  cylindrical  coordinate  system 
for  both  vector(momentum)  field  and  scalar  (temperature,  concentration,  etc.)  field  by 
introducing  source-like  term  and  body  force-like  term(6)  in  lattice  Boltzmann  equations. 
The  model  has  the  second  order  accuracy  in  space  and  is  easy  to  be  implemented. 

The  remainder  of  this  Appendix  is  organized  as  follows.  Lattice  Boltzmann  equations 
for  axisymmetric  flows  in  cylindrical  coordinate  system  are  given  in  Section  G.2,  in  which 
macroscopic  equations  such  as  continuity  equation,  momentum  equations,  scalar  equation 
for  temperature  or  concentration  are  derived;  Numerical  evaluation  are  presented  in  Section 
G.3;  Finally,  we  conclude  in  Section  G.4  with  a  brief  discussion. 

G.2  Lattice  Boltzmann  model  for  axisymmetric  flows 

G.2.1  Hydrodynamic  equations  in  cylindrical  coordinate  system 

In  cylindrical  coordinate,  continuity  equation,  Navier-stokes  equations  and  scalar(temperature, 
mass  concentration,  etc.)  equation  for  incompressible  flow  are  expressed  as  following 
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dp  1  djprur)  1  djpug)  d(puz) 

dt^  r  dr  +  r  dO  +  ~dT~  ~  VG-1) 


dur  duT  no  duT 
dt  ~irUr  dr  +  r  dO 


dur 

lh' 


T 


1  dp  ,  iA&Ut  ,  1  dur  ,  1  d2ur  t  d2ur  uT  |  2  dug  1 
pdr  1  dr2  r  dr  r2  d62  +  dz 2  r2+r2  dO  * 

(G.2) 


dug  dug  Ug  dug  dug  UgUr 

-5-+Ur— H — -+u  — _+ 

dt  dr  r  80  dz  r 


1  |  1  ^9  1  Pug  ue2dur 

pr  d6  dr 2  r  dr  r2  d2  dz2  r2"*"r2  dd 

(G.3) 


duz  duz  ug  duz  duz 
dt  +Ur  dr  +T~de+Uz~d7 


,  1  ^u,  , 

pdz  1  dr2  r  dr  r2  d02  +  dz2  J 


(G.4) 


dF  dF  ug  dF  dF 

Q.  '  r\  +  o/3  'U'Z’Tl 

dt  dr  r  dO  dz 


.  A  d  .  dF.  1  d2F  d2F .  „ 

+  ^~dW  +  dZ^]  +  Qp 


(G.5) 


where  v  is  the  kinetic  viscosity  of  the  fluid,  Ap  is  the  scalar  related  property  such  as  the 
thermal  drffusivity  for  temperature  and  diffusivity  coefficient  for  concentration,  and  Qp  is 
a  source  term  for  the  scalar  field. 


In  an  axisymmetric  field  with  ug  =  0,  above  equations  can  be  simplified  as 

dp  dpur  dpuz  pur 

dt+  dr  +  dz  ~  (°-6) 


dur  1  dp  o  V  dur  ur  N 


(0.7) 


and 


duz  1  dp  2  v  duz 

—  +  n^u,~-~  +  ^u,  +  -w. 


dF 

~  +  u-  VF  =  AfV2F  +  sF. 


(G.8) 


(G.9) 


respectively  with  u  =  ureT  +  uzez,V  =  dTer  +  8zez  ,  V2  =  d2  +  d2,  and  sF  =  +  QF. 
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G.2.2  Lattice  Boltzmann  equation  for  vector  (momentum)  field 


Comparing  Eqs.  (G.6)  and  (G.8)  with  2D  continuity  equation  and  momentum  equations 
in  Cartesian  coordinate  respectively,  we  found  there  is  an  extra  term  of  s  =  —^r  in  Eq. 
(G.6)  and  a*  =  ~{drUi  —  i  =  r, z  in  Eq.  (G.8).  This  suggests  us  to  introduce  a  term 
ga(r,t)  in  lattice  Boltzmann  equation(7;  8) 


/■»(>■  +  ea5t,  t  +  S,)  =  fjr,  t)  --[fa-  f«q’}  +  Sisa(r,  i), 


(G.10) 


where  fa  is  the  density  distribution  function  with  discrete  velocity  ea  along  ath  direction, 
fa^  is  the  equilibrium  distribution  function,  and  r  is  the  relaxation  time  determining  the 
viscosity  v  of  the  modeled  fluid.  We  use  the  nine  velocity  model  in  two  dimensions  (D2Q9) 
in  our  analysis.  The  equilibria  for  incompressible  flow(10)  are 


ft* = {/>' 


+  Po 


3ea  •  u  9(eQ  ■  u )2 


3u2 

2c2 


}■ 


c‘  2c4 

where  p 1  is  the  density  fluctuation,  and  po  is  the  mean  density  in  the  system  which  is  usually 
set  to  one.  The  total  density  is  p  =  pQ  +  p' .  The  mass  and  momentum  conservations  are 
enforced: 


»'=£*>  =  Eft*’ 

a  a 

p0U  =  eafa  =  ea/ieq)- 
a  Cl 

The  spatial  and  velocity  dependent  microscopic  term  ga(r ,  t)  is  to  reflect  two  extra  terms 
mentioned  above,  source-like  term  s  in  Eq.  (G.6)  and  body-force-like  term  a  in  Eq.  (G.8) 
simultaneously.  We  shall  use  Chapman- Enskog  technique  to  determine  it. 

Chapman-Enskog  expansions  of  fa  and  dt  in  the  powers  of  e  are  the  same  as  Eq.(C.16). 
For  purpose  of  extracting  the  dynamics  of  this  model,  we  perform  a  chapman-Enskog  type 
expansion  with  the  ga  as  well  as  fa, 


(G.U) 
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(G.12) 


Here  we  have  assumed  that  there  is  no  zeroth  term  in  the  expansion(ll). 

Taylor  expansion  in  St  to  the  second  order  and  a  substitution  of  chapman-Enskog  ex¬ 


pansions  lead  to 


where 


fa(x  +  ea5t,t  +  St)  —  fa{x,t) 

=  ftDtU  +  +  o((Je)3) 

=  eS,D,jM  +  ^(9,,/f  +  +  o(e3), 


Dt  =  dt  +  ea-  V,  Dt0  =  dt0  +  eQ  •  V. 


The  first  few  equations  of  a  set  of  successive  equations  in  the  order  of  e  obtained  from  Eq. 
(G.10)  are 

=/<*»,  (G.13) 


*‘-A  ,fSll  =  ~fP+s2\ 


r2  .  y(0)  _|_  (2t  1  ]Dtofa  ^  _  StDtogi  )  —  fg  )  ^  ^(2) 


(G.13) 

(G.14) 


(G.15) 


The  distribution  function  fa  is  the  normal  solution  which  is  constrained  by 


£/f 


£/in) 


n  >  0 


(G.16) 


(G.17) 


For  D2Q9  model,  the  tensor 


g(n)  _  'Sp  u;aeQ)iea)2  •  •  •  ec 


(G.18) 
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dtj  +  V  ■  (p0u)  =  s  (G.30) 

dtoiPoUi )  +  +  pou  •  Vm  =  0  (G.31) 

and 

dhP'  =  E^2)  ~\&tdtoS  (G.32) 

a 

2r _ ^  2 

dti  (poui) - g — dtc2po(dl  +  dy  )ui  =  ^  eai9^  —  (G.33) 

a 

Combining  the  first  and  the  second  order  of  e  above,  we  obtain 

dtp'  +  V  •  (p0«)  =  3  +  -  ^5tdtos,  (G.34) 

a 

+  u  •  Vui  +  -  i/(fl“  +  2  e«i^2)  -  ^—StrdiS,  (G.35) 

Po  po  a  Opo 

where  p  =  c2p’  and  v  =  (2r  -  1  )c25tpo/6. 

For  the  purpose  of  recovering  the  macroscopic  equations  (G.6)  and  (G.8),  we  have  the 
following  constrains  for  p£2)  by  comparing  Eq.  (G.34)  with  Eq.  (G.6)  and  Eq.(G.35)  with 
Eq.  (G.8)  on  the  right-hand  side  respectively, 

J^9(a}  =  \stStos  =  s',  (G.36) 

OC 

£»2  2 

^  eod9a^  =  POai  +  -r^tT^iS  =  P°ai  +  Ta'v 
a  6  6 

where  a[  =  dts.  These  equalities  are  achieved  if  we  select 

9a  *  =  £(was)  +  £2(^j-waea  ■  a  +  waea  ■  a')  (G.37) 


Notice  from  Eqs.  (G.32)  and  (G.33)  that 
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Finally,  ga  can  be  expressed  as 


_  .  3po  1  c2  dp'  po  . 

9a  —  was  H — 2"iyaeQ  •  a  +  -- — —  +  — u  •  Vur  +  wQea  •  a 
c  l  or  or  r 


(G.38) 


with 


s  =  — 


PoV>r 


uduz  v  ,dur  ur. 

az  =  — X-,  ar  =  -(— -  -  -L , 
r  or  r  dr  r 


,  ds  ,  ds 

a*  =  StrTz’  ar  =  s,Tfr 

Substitute  Eq.  (G.37)  into  Eqs.  (G.34)  and  (G.35),  the  continuity  equation  (G.6)  and 
incompressible  Navier-Stokes  equation  (G.8)are  recovered. 


G.2.3  Lattice  Boltzmann  equation  for  scalar  field 

Following  the  same  methodology  as  LBM  to  vector  field,  the  lattice  Boltzmann  equation 
for  scalar  F  is  written  as 


Fa{r  +  eaSt,t  +  St)  =  Fa(r,t) - [FQ  —  F^]  +  Stsa(r,t), 


(G.39) 


where  sa  =  waSF 

The  equilibria  for  D2Q9  mesh  are 


4eq)  =  waF[  1  +  +  'f? 

0  LC 


3u2 

2?J 


Conservation  law  corresponding  to  the  scalar  considered  is  enforced 


(0.40) 

a  a 

Taylor  expanding  to  the  second  order  of  St  yields 

Fair  +  +  St)  —  Fa(r,t)  =  StDfFa  +  - StDtDfFa  4 -  (G.41) 

Scalar  transport  (diffusion,  thermal  conduction,  etc.)  is  a  slow  process  on  large  spatial 
scales  (12),  which  suggests  the  following  Chapman-Enskog  expansion 
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Fa  =  Fi0)  +  eF^  +  £2fW  +  o(e3)  (G.42) 

dt  =  e2dh  +  o(e3)  (G.43) 

Substitutions  of  Taylor  expansion  (G.41)  and  Chapman-Enskog  expansions  (G.42)  and 
(G.43)  into  Eq.  (G.39)  bring 

e(e  •  V)4°>  +  e2[dhF( <°>  +  (i  -  r)St(e  ■  V)2fW]  +  o(e 3) 

=  +  e*F, <2)  +  o(e3))  +  sa 

Summing  over  all  a  of  above  equation,  Eq.  (G.9)  is  recovered  by  noticing  fP  =  Fp^  (so 

that  E  FP  =  E  Fi2)  =  0)  and  fP  =  -r8t(e  •  V)Fi0)  with  A.p  = 

a  a  b 


G.3  Evaluation  simulations 


In  order  to  evaluate  the  lattice  Boltzmann  models  presented  above,  we  perform  simulations 
of  two  problems  as  following. 

To  discrete  V xa  and  V p  in  the  binary  interaction  term  in  Eq.  (C.9)  and  Vur  in  Eq. 
(G.38),  we  use  the  following  stencil  for  dXl  and  dX2: 


-1 

0 

1 

«  1 
’  ^“12 

1 

4 

1 

-4 

0 

4 

0 

0 

0 

1 - 

1 

0 

1 

-1 

-4 

- 1 

i — I 

1 

dXlf(x)  »  i[4/(x  +  e1)  +  /(f+e5)  +  /(£+e8) 

-  f(x  +  e6)-4f(x  +  e3)-f(x  +  e7)] 

9x2 /(®)  ~  n[4/(x  +  e-j)  +  /(£  +  eg)  +  f(x  +  e§) 

-  f(x  +  e7)-4f(x  +  e4)-f(x  +  e&)} 


(G.44) 


(G.45) 


(G.46) 


in  which  /  can  be  xc  and  pa  or  ur  and  uz. 
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G.3.1  Axisymmetric  flows 


In  this  part,  our  exercises  axe  pipe  flow  and  counter  flow  by  using  the  model  for  vector  field. 
Fully  developed  pipe  flow:  In  a  pipe  with  the  radius  R,  flow  is  driven  by  a  constant 
pressure  gradient  G  =  —dp/dz.  The  geometric  configuration  of  the  pipe  flow  is  shown  in 
Fig.  (G.l). 

Assuming  that  the  flow  is  a  fully  developed  laminar  flow,  the  analytical  solution  of 
velocity  profile  across  the  flow  stream  is 

Uz{r)  =  ~^~o{r2  ~  R2)-  (G‘47) 

Our  simulation  is  carried  out  in  the  shadowed  area  of  Fig.  (G.l).  Non-slip  and  symmetric 
boundary  conditions  are  applied  at  the  up  and  down  parallel  walls  respectively  and  periodic 
boundary  at  inlet  and  outlet.  The  constant  pressure  gradient  is  treated  as  body  force(13). 
The  system  size  is  Nz  x  Nr  =  250  x  50.  The  initial  condition  is  set  as  u  =  0,  po  =  1,  and 
p'  =  0  everywhere  in  the  domain.  Two  cases  axe  run  with  G  =  2.89351e  —  05, r  =  1.0  and 
G  =  1.73611e  —  06, r  =  0.8  respectively. 

The  velocity  profile  comparisons  between  analytical  solution  from  Eq.  (G.47)  and 
LBM  simulation  solution  are  shown  in  Fig.  (G.2).  Figures  show  that  simulation  solu- 
tions(symbols)  axe  in  excellent  agreement  with  analytical  solutions. 

The  L2  norm  error  between  analytical  solution  and  LBM  simulation  solution  is  defined 


as 


_E2  = 


(G.48) 


In  figure  (G.3),  the  relation  of  E2  with  the  grid  size  h  =  l/(Nr  —  2)  shows  that  current 
model  has  second  order  accuracy  in  space. 


Counterflow:  Counterflow  field  is  established  by  impinging  two  symmetrical  counterflow¬ 
ing  streams  from  two  identical  nozzles  facing  up  and  down(see  Fig.(G.4)). 
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The  simulation  is  in  shadowed  area.  Initially,  we  set  po  ,  p',and  u  equal  to  zero  every¬ 
where.  Symmetric  boundary  conditions  are  applied  at  left  and  bottom  boundaries  and  fully 
developed  boundary  at  outlet.  A  uniform  velocity  is  given  at  inlet  (z  =  —  Uin  and  ur  =  0). 
Grid  resolution  is:  Nr  x  Nz  =  300  x  180.  We  use  r  =  1.0,  U{n  =  0.1  in  simulation. 

Figure  (G.5)  shows  velocity  field  after  steady  state  is  achieved.  In  Figure  (G.6),  M  is 
the  instant  total  mass  while  Mq  is  the  initial  total  mass  in  the  simulation  domain.  Figure 
(G.6)  shows  the  mass  conservation  in  the  domain  after  a  steady  flow  is  achieved. 

G.4  Conclusions 

We  have  presented  a  LBM  model  for  axisymmetric  flows  in  cylindrical  coordinate  system. 
Macroscopic  equations  governing  the  flows  are  derived  in  details  by  the  Chapman-Enskog 
technique. 

Simulations  using  this  model  are  optimistic.  The  scheme  has  been  proved  in  the  sec¬ 
ond  order  of  accuracy.  In  the  fully  developed  pipe  flow  simulation,  velocity  profiles  were 
in  excellent  agree  with  analytical  predictions;  while  in  the  counterflow  simulation,  mass 
conservation  was  verified  and  stable  velocity  field  was  obtained. 
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G.5  Figures 


Figure  G.l:  A  schematic  illustration  of  pipe  flow  driven  by  a  constant  pressure  gradient  G. 
Only  the  shadowed  region  is  simulated  according  to  the  symmetry  feature. 


Figure  G.2:  Velocity  profiles  of  a  fully  developed  pipe  flow.  Lines:  analytical  solutions; 
Symbols:  LBM  simulation  solutions,  (a)  G  =  2.89351e  —  05, r  =  1.0;  (b)  G  =  1.73611e  - 
06, t  =  0.8 
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